1 Task info and setup

This R Markdown script analyses data from the PAL (probabilistic associative learning) task of the EMBA project. It focuses on the comparison of adults with ADHD (with or without comorbid ASD) and adults without psychiatric diagnoses. The data was preprocessed before being read into this script.

In this task, participants view faces for which they have to judge whether the portrayed emotion is positive or negative. The faces are preceded by a tone which is predictive of the valence of the following face. The visual angle of the faces was 3.46 degrees wide and 4.72 degrees high.

1.1 Plot the development of associations across task

# load the trial order
df.trl = read_csv('../data/PAL_scheme.csv', show_col_types = F) %>%
  mutate(
    association = ut*100, 
    prob_HP = prob_HP*100
  )

# create rectangles to colour the phases: prevolatile, volatile and postvolatile
rects = data.frame(xstart = c(-Inf, 72, 264), 
                   xend = c(72, 264, Inf), 
                   col = c("a", "b", "c"))
  
ggplot() +
  geom_rect(data = rects, aes(xmin = xstart, xmax = xend, 
                              ymin = -Inf, ymax = Inf, 
                              fill = col), alpha = 0.4) +
  scale_fill_manual(values=c(c_mid_highlight, c_green, c_dark_highlight), 
                    guide = "none") +
  geom_jitter(data = df.trl, 
              aes(x = trl, y = association, shape = expected),
              colour = "black", alpha = 0.8, size = 1, width = 0, height = 4) +
  geom_line(data = df.trl, 
            aes(x = trl, y = prob_HP),
            colour = "black") +
  scale_shape_manual(values = c(8, 1, 3)) +
  theme_bw() + 
  scale_y_continuous(breaks=c(16.7, 83.3)) +
  labs(x = "trial", y = "percent", 
       title = "High tone > positive, low tone > negative emotion") +
  annotate("text", x = 36, y = 125, label= "prevolatile", size = 4) +
  annotate("text", x = 168, y = 125, label= "volatile", size = 4) +
  annotate("text", x = 301, y = 125, label= "postvolatile", size = 4) +
  theme(legend.position = c(0.9, 0.5), legend.title = element_blank(), 
        legend.background = element_rect(fill = alpha("white", 0)),
        legend.key = element_rect(fill = alpha("white", 0)),
        plot.title = element_text(hjust = 0.5), text = element_text(size = 12))

ggsave("plots/PAL_scheme.png", 
       units = "cm",
       width = 27,
       height = 9)

1.2 Some general settings

# number of simulations
nsim = 250

# set number of iterations and warmup for models
iter = 6000
warm = 1500

# set the seed
set.seed(2468)

1.3 Package versions

The following packages are used in this RMarkdown file:

## [1] "R version 4.5.1 (2025-06-13)"
## [1] "knitr version 1.50"
## [1] "ggplot2 version 3.5.2"
## [1] "brms version 2.22.0"
## [1] "bridgesampling version 1.1.2"
## [1] "tidyverse version 2.0.0"
## [1] "ggpubr version 0.6.1"
## [1] "ggrain version 0.0.4"
## [1] "bayesplot version 1.13.0"
## [1] "SBC version 0.3.0.9000"
## [1] "rstatix version 0.7.2"
## [1] "officer version 0.6.10"
## [1] "effectsize version 1.0.1"
## [1] "BayesFactor version 0.9.12.4.7"
## [1] "bayestestR version 0.17.0"

1.4 Analysis guidelines

We planned to determine the group-level effect subjects following Barr (2013). For each model, experiment specific priors were set based on previous literature or the task (see comments in the code).

We performed prior predictive checks as proposed in Schad, Betancourt and Vasishth (2020) using the SBC package before running this analysis (see [!LINK]). We performed these data-naive evaluations based on the comparison of three groups. To do so, we create 250 simulated datasets where parameters are simulated from the priors. These parameters are used to create one fake dataset. Both the true underlying parameters and the simulated discrimination values are saved. Then, we create graphs showing the prior predictive distribution of the simulated discrimination threshold to check whether our priors fit our general expectations about the data. Next, we perform checks of computational faithfulness and model sensitivity as proposed by Schad, Betancourt and Vasishth (2020) and implemented in the SBC package. We create models for each of the simulated datasets. Last, we calculate performance metrics for each of these models, focusing on the population-level parameters.

1.5 Preparation

First, we load the data and combine it with demographic information including the diagnostic status of the subjects. Then, all predictors are set to sum contrasts.

# check if the data file exists, if yes load it:
if (!file.exists("../data/PAL-ADHD_data.RData")) {
  
  path = "/home/emba/Documents/EMBA"

  # get demo info for subjects
  df.sub = read_csv(file.path(path, "CentraXX", "EMBA_centraXX.csv"), 
                    show_col_types = F) %>%
    mutate(
      diagnosis = as.factor(recode(diagnosis, "CTR" = "COMP"))
    )
  
  # set the data path
  dt.path  = file.path(path, "BVET")
  dt.explo = file.path(path, "BVET-explo")
  
  # load the preprocessed data: df.tsk
  df.tsk = rbind(readRDS(file.path(dt.path,  "PAL_tsk.rds")),
                 readRDS(file.path(dt.explo, "PAL_tsk.rds")))
  
  # load excluded participants (accuracy < 2/3)
  exc = c(scan(file.path(dt.path, 'PAL_exc.txt'),  what="character", sep=NULL),
          scan(file.path(dt.explo, 'PAL_exc.txt'), what="character", sep=NULL))
  
  # merge with group and focus on comparison and autistic adults
  df.tsk = merge(df.sub %>% select(subID, diagnosis), 
                 df.tsk, all.y = T) %>%
    filter(diagnosis != "ASD" & !(subID %in% exc)) %>%
    droplevels() %>%
    mutate_if(is.character, as.factor)
  
  # only keep participants included in the study in the subject data frame
  df.sub = merge(df.tsk %>% select(subID) %>% distinct(), df.sub, all.x = T) %>%
    droplevels()

  # anonymise the data based on predetermined file > same for all data
  df.recode = read_csv(file.path(path, "BVET", "PID_anonymisation.csv"))
  recode = as.character(df.recode$subID)
  names(recode) = df.recode$PID
  df.tsk$subID = str_replace_all(df.tsk$subID, recode)
  
  # save the data for use with the EMBA_HGF toolbox
  df.tsk %>% select(subID, trl, diagnosis, ut, difficulty, acc, rt) %>%
    ungroup() %>%
    arrange(subID, trl) %>%
    write_csv(., file = "../data/PAL-ADHD_data.csv")
  
  # print gender frequencies and compare them across groups
  tb.gen = xtabs(~ gender + diagnosis, data = df.sub)
  ct.full = contingencyTableBF(tb.gen, 
                               sampleType = "indepMulti", 
                               fixedMargin = "cols")
  
  # get all self-descriptions for DAN group
  gen.desc = unique(tolower(df.sub[df.sub$gender == "dan",]$gender_desc))
  
  # create a demographics overview
  ls.demo = createDFdemo(df.sub, 
                         c("age", "edu", "BDI_total", "ASRS_total", "RAADS_total", "TAS_total", "iq"),
                         "diagnosis")
  
  # save it all
  save(df.tsk, ls.demo, ct.full, gen.desc, tb.gen, 
       file = "../data/PAL-ADHD_data.RData")
  
} else {
  
  load("../data/PAL-ADHD_data.RData")
  
}

# print the outcome of the two contingency tables for comparison: all participants
ct.full@bayesFactor
##                         bf error                     time         code
## Non-indep. (a=1) -2.156359     0 Fri Sep  5 12:54:03 2025 21ea4a28ddd8
# print gender and diagnosis distribution
kable(tb.gen)
ADHD BOTH COMP
dan 2 3 0
fem 9 12 10
mal 11 7 12
# drop the neutral condition for the analysis
df.pal = df.tsk %>%
  filter(expected != "neutral" & !is.na(rt.cor)) %>% droplevels() %>%
  mutate_if(is.character, as.factor) %>%
  ungroup()

# calculate variance of reaction times: 
# no difficulty due to suboptimal posterior fit
df.var = df.tsk %>%
  filter(expected != "neutral") %>% droplevels() %>%
  group_by(subID, diagnosis, phase, expected) %>% 
  summarise(
    totaln = n(),
    valuen = sum(!is.na(rt.cor)),
    rt.var = sd(rt.cor, na.rm = T)
  ) %>%
  mutate(
    perc = valuen/totaln
  ) %>% filter(perc >= 2/3) %>%
  mutate_if(is.character, as.factor) %>% ungroup()

# set and print the contrasts
contrasts(df.pal$diagnosis) = contr.sum(3)
contrasts(df.pal$diagnosis)
##      [,1] [,2]
## ADHD    1    0
## BOTH    0    1
## COMP   -1   -1
contrasts(df.pal$expected) = contr.sum(2)
contrasts(df.pal$expected)
##            [,1]
## expected      1
## unexpected   -1
contrasts(df.pal$phase) = contr.sum(3)
contrasts(df.pal$phase)
##              [,1] [,2]
## prevolatile     1    0
## volatile        0    1
## postvolatile   -1   -1
contrasts(df.pal$difficulty) = contr.sum(3)
contrasts(df.pal$difficulty)
##           [,1] [,2]
## easy         1    0
## medium       0    1
## difficult   -1   -1
contrasts(df.var$diagnosis) = contr.sum(3)
contrasts(df.var$diagnosis)
##      [,1] [,2]
## ADHD    1    0
## BOTH    0    1
## COMP   -1   -1
contrasts(df.var$expected) = contr.sum(2)
contrasts(df.var$expected)
##            [,1]
## expected      1
## unexpected   -1
contrasts(df.var$phase) = contr.sum(3)
contrasts(df.var$phase)
##              [,1] [,2]
## prevolatile     1    0
## volatile        0    1
## postvolatile   -1   -1
# print final group comparisons for the paper into a Word document
read_docx() %>%
    body_add_table(ls.demo$df.demo %>% arrange(measurement) %>% 
                         mutate(bf.log = 
                                  if_else(
                                    bf.log > log(3), 
                                    sprintf("%.3f*", bf.log),
                                    sprintf("%.3f", bf.log)))) %>%
    print(target = "PAL-ADHD_demo.docx")

read_docx() %>%
    body_add_table(ls.demo$df.post %>% arrange(measurement) %>% 
                         mutate(bf.log = 
                                  if_else(
                                    bf.log > log(3), 
                                    sprintf("%.3f*", bf.log),
                                    sprintf("%.3f", bf.log)))) %>%
    print(target = "PAL-ADHD_post.docx")

2 Reaction time variance

In the preregistration, we noted the following population-level effects for the model of the reaction time variances: group, expectancy, phase and difficulty. However, the posterior fit for this model was suboptimal, therefore, we dropped the predictor difficulty.

2.1 Model setup

# figure out slopes for subjects
kable(head(df.var %>% count(subID, expected)))
subID expected n
1 expected 3
1 unexpected 3
10 expected 3
10 unexpected 3
11 expected 3
11 unexpected 3
kable(head(df.var %>% count(subID, phase)))
subID phase n
1 prevolatile 2
1 volatile 2
1 postvolatile 2
10 prevolatile 2
10 volatile 2
10 postvolatile 2
kable(head(df.var %>% count(subID, expected, phase)))
subID expected phase n
1 expected prevolatile 1
1 expected volatile 1
1 expected postvolatile 1
1 unexpected prevolatile 1
1 unexpected volatile 1
1 unexpected postvolatile 1
# model formula
f.var = brms::bf(rt.var ~ diagnosis * expected * phase +
                   (expected + phase | subID)
                 )

# set weakly informative priors
priors = c(
  prior(normal(4.5,  0.50),  class = Intercept),
  prior(normal(0,    0.50),  class = sigma),
  prior(normal(0.15, 0.15),  class = sd),
  prior(normal(0,    0.25),  class = b),
  prior(lkj(2),              class = cor)
)

2.2 Posterior predictive checks

As the next step, we fit the model, check whether there are divergence or rhat issues, and then check whether the chains have converged.

# fit the final model
m.var = brm(f.var, seed = 4646,
            df.var, prior = priors,
            family = lognormal,
            iter = iter, warmup = warm,
            backend = "cmdstanr", threads = threading(8),
            file = file.path(brms_dir, "m_var"),
            save_pars = save_pars(all = TRUE)
            )
rstan::check_hmc_diagnostics(m.var$fit)
## 
## Divergences:
## 0 of 18000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 18000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.
# check that rhats are below 1.01
sum(brms::rhat(m.var) >= 1.01, na.rm = T)
## [1] 0
# check the trace plots
post.draws = as_draws_df(m.var)
mcmc_trace(post.draws, regex_pars = "^b_",
           facet_args = list(ncol = 4)) +
  scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
  scale_y_continuous(breaks=scales::pretty_breaks(n = 3))
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.

This model has no pathological behaviour with E-BFMI, no divergent sample and no rhats that are higher or equal to 1.01. Therefore, we go ahead and perform our posterior predictive checks.

# get posterior predictions
post.pred = posterior_predict(m.var, ndraws = nsim)

# check the fit of the predicted data compared to the real data
p1 = pp_check(m.var, ndraws = nsim) + 
  theme_bw() + theme(legend.position = "none") + xlim(0, 1000)

# distributions of means compared to the real values per group
p2 = ppc_stat_grouped(df.var$rt.var, post.pred, df.var$diagnosis) + 
  theme_bw() + theme(legend.position = "none")

# distributions of means compared to the real values per phase
p3 = ppc_stat_grouped(df.var$rt.var, post.pred, df.var$phase) + 
  theme_bw() + theme(legend.position = "none")

# distributions of means compared to the real values per expected
p4 = ppc_stat_grouped(df.var$rt.var, post.pred, df.var$expected) + 
  theme_bw() + theme(legend.position = "none")

p = ggarrange(p1, p2, p3, p4,
          nrow = 4, ncol = 1, labels = "AUTO")
annotate_figure(p, top = text_grob("Posterior predictive checks: RT variance", 
                                   face = "bold", size = 14))

The simulated data based on the model (light blue) fits well with the real data (dark blue), both for the overall distribution (A) and the aggregated scores for the diagnostic groups (B).

2.3 Inferences

Now that we are convinced that we can trust our model, we have a look at its estimate and use the hypothesis function to assess our hypotheses and perform explorative tests.

# print a summary
summary(m.var)
##  Family: lognormal 
##   Links: mu = identity; sigma = identity 
## Formula: rt.var ~ diagnosis * expected * phase + (expected + phase | subID) 
##    Data: df.var (Number of observations: 387) 
##   Draws: 4 chains, each with iter = 6000; warmup = 1500; thin = 1;
##          total post-warmup draws = 18000
## 
## Multilevel Hyperparameters:
## ~subID (Number of levels: 66) 
##                          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept)                0.36      0.03     0.30     0.43 1.00     3864
## sd(expected1)                0.02      0.01     0.00     0.05 1.00     4508
## sd(phase1)                   0.09      0.02     0.05     0.12 1.00     6414
## sd(phase2)                   0.03      0.02     0.00     0.08 1.00     4058
## cor(Intercept,expected1)    -0.20      0.31    -0.74     0.47 1.00    22631
## cor(Intercept,phase1)       -0.23      0.18    -0.56     0.13 1.00    17311
## cor(expected1,phase1)        0.26      0.33    -0.48     0.81 1.00     2930
## cor(Intercept,phase2)        0.04      0.30    -0.56     0.62 1.00    22936
## cor(expected1,phase2)        0.19      0.37    -0.59     0.81 1.00     6950
## cor(phase1,phase2)           0.08      0.34    -0.58     0.70 1.00    17998
##                          Tail_ESS
## sd(Intercept)                7484
## sd(expected1)                5809
## sd(phase1)                   6234
## sd(phase2)                   6446
## cor(Intercept,expected1)    12093
## cor(Intercept,phase1)       13589
## cor(expected1,phase1)        5561
## cor(Intercept,phase2)       12392
## cor(expected1,phase2)       11325
## cor(phase1,phase2)          14323
## 
## Regression Coefficients:
##                             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept                       4.71      0.04     4.62     4.79 1.00     2209
## diagnosis1                      0.09      0.06    -0.03     0.21 1.00     2291
## diagnosis2                      0.03      0.06    -0.09     0.15 1.00     2444
## expected1                       0.06      0.01     0.04     0.08 1.00    30294
## phase1                          0.02      0.02    -0.01     0.05 1.00    15167
## phase2                          0.03      0.01     0.00     0.06 1.00    26066
## diagnosis1:expected1           -0.01      0.01    -0.04     0.02 1.00    22496
## diagnosis2:expected1           -0.00      0.01    -0.03     0.02 1.00    22150
## diagnosis1:phase1               0.00      0.02    -0.05     0.05 1.00    13212
## diagnosis2:phase1               0.05      0.02    -0.00     0.10 1.00    13450
## diagnosis1:phase2              -0.03      0.02    -0.07     0.01 1.00    19605
## diagnosis2:phase2              -0.02      0.02    -0.06     0.02 1.00    20126
## expected1:phase1                0.02      0.01    -0.00     0.05 1.00    27027
## expected1:phase2               -0.02      0.01    -0.05     0.00 1.00    25645
## diagnosis1:expected1:phase1    -0.01      0.02    -0.05     0.03 1.00    19059
## diagnosis2:expected1:phase1    -0.01      0.02    -0.05     0.02 1.00    19306
## diagnosis1:expected1:phase2     0.04      0.02    -0.00     0.07 1.00    19797
## diagnosis2:expected1:phase2     0.00      0.02    -0.03     0.04 1.00    18557
##                             Tail_ESS
## Intercept                       5022
## diagnosis1                      4509
## diagnosis2                      4865
## expected1                      13918
## phase1                         13959
## phase2                         13147
## diagnosis1:expected1           14370
## diagnosis2:expected1           14145
## diagnosis1:phase1              13651
## diagnosis2:phase1              13448
## diagnosis1:phase2              13824
## diagnosis2:phase2              14268
## expected1:phase1               15127
## expected1:phase2               15707
## diagnosis1:expected1:phase1    15059
## diagnosis2:expected1:phase1    13226
## diagnosis1:expected1:phase2    14941
## diagnosis2:expected1:phase2    14865
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.18      0.01     0.17     0.20 1.00     6545    10033
## 
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# get the estimates and compute group comparisons
df.m.var = post.draws %>% 
  select(starts_with("b_")) %>%
  mutate(
    b_COMP         = - b_diagnosis1 - b_diagnosis2,
    b_postvolatile = - b_phase1 - b_phase2,
    BOTH           = b_Intercept + b_diagnosis2,
    ADHD           = b_Intercept + b_diagnosis1,
    COMP           = b_Intercept + b_COMP,
    `ADHD-BOTH`    = ADHD - BOTH,
    `ADHD-COMP`    = ADHD - COMP,
    `BOTH-COMP`    = BOTH - COMP
  )

# plot the posterior distributions
df.m.var %>%
  select(starts_with("b_")) %>%
  pivot_longer(cols = starts_with("b_"), names_to = "coef", values_to = "estimate") %>%
  subset(!startsWith(coef, "b_Int")) %>%
  mutate(
    coef = substr(coef, 3, nchar(coef)),
    coef = str_replace_all(coef, ":", " x "),
    coef = str_replace_all(coef, "diagnosis1", "ADHD"),
    coef = str_replace_all(coef, "diagnosis2", "ADHD+ASD"),
    coef = str_replace_all(coef, "expected1", "expected"),
    coef = str_replace_all(coef, "expected2", "unexpected"),
    coef = str_replace_all(coef, "phase1", "prevolatile"),
    coef = str_replace_all(coef, "phase2", "volatile"),
    coef = fct_reorder(coef, desc(estimate))
  ) %>% 
  group_by(coef) %>%
  mutate(
    cred = case_when(
      (mean(estimate) < 0 & quantile(estimate, probs = 0.975) < 0) |
        (mean(estimate) > 0 & quantile(estimate, probs = 0.025) > 0) ~ "credible",
      T ~ "not credible"
    )
  ) %>% ungroup() %>%
  ggplot(aes(x = estimate, y = coef, fill = cred)) +
  geom_vline(xintercept = 0, linetype = 'dashed') +
  ggdist::stat_halfeye(alpha = 0.7) + ylab(NULL) + theme_bw() +
  scale_fill_manual(values = c(c_dark, c_light)) + theme(legend.position = "none")

# H1a: COMP < ADHD
h1a = hypothesis(m.var, "0 < 2*diagnosis1 + diagnosis2")
h1a
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (0)-(2*diagnosis1... < 0    -0.21      0.11    -0.39    -0.03      37.54
##   Post.Prob Star
## 1      0.97    *
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
## Exploration

# COMP < BOTH
e1 = hypothesis(m.var, "0 < diagnosis1 + 2*diagnosis2", alpha = 0.025)
e1
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (0)-(diagnosis1+2... < 0    -0.15      0.11    -0.37     0.06      11.82
##   Post.Prob Star
## 1      0.92     
## ---
## 'CI': 95%-CI for one-sided and 97.5%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 97.5%;
## for two-sided hypotheses, the value tested against lies outside the 97.5%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
# ADHD vs. BOTH
e2 = hypothesis(m.var, "diagnosis1 > diagnosis2", alpha = 0.025)
e2
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (diagnosis1)-(dia... > 0     0.06      0.11    -0.15     0.27       2.53
##   Post.Prob Star
## 1      0.72     
## ---
## 'CI': 95%-CI for one-sided and 97.5%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 97.5%;
## for two-sided hypotheses, the value tested against lies outside the 97.5%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
## Exploration of possible transition effects

## create design matrix to figure out how to set comparisons for hypotheses
df.des = cbind(df.var, 
               model.matrix(~ diagnosis * expected * phase, 
                            data = df.var)) %>%
  select(-subID, -rt.var, -perc, -totaln, -valuen) %>% distinct()

# COMP(vol-pre) > ADHD(vol-pre)
as.data.frame(t(df.des %>% 
    ungroup() %>%
    filter((diagnosis == "ADHD" | diagnosis == "COMP") & 
             ((phase == "volatile") | (phase == "prevolatile"))) %>%
    group_by(diagnosis, phase) %>%
    summarise(across(where(is.numeric), ~ mean(.x))) %>%
    group_by(diagnosis) %>%
    summarise(across(where(is.numeric), ~ diff(.x))) %>%
    select(where(is.numeric)) %>%
    map_df(~ diff(.x)))) %>%
  filter(V1 != 0)
##                   V1
## diagnosis1:phase1  2
## diagnosis2:phase1  1
## diagnosis1:phase2 -2
## diagnosis2:phase2 -1
e3 = hypothesis(m.var, "0 < 2*diagnosis1:phase1 + diagnosis2:phase1 + 
                - (2*diagnosis1:phase2 + diagnosis2:phase2)", 
                alpha = 0.025)
e3
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (0)-(2*diagnosis1... < 0    -0.13      0.06    -0.26    -0.01      46.62
##   Post.Prob Star
## 1      0.98    *
## ---
## 'CI': 95%-CI for one-sided and 97.5%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 97.5%;
## for two-sided hypotheses, the value tested against lies outside the 97.5%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
# COMP(vol-pre) > BOTH(vol-pre)
e4 = hypothesis(m.var, "0 < diagnosis1:phase1 + 2*diagnosis2:phase1 +
                - (diagnosis1:phase2 + 2*diagnosis2:phase2)", 
                alpha = 0.025)
e4
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (0)-(diagnosis1:p... < 0    -0.17      0.06    -0.29    -0.04     201.25
##   Post.Prob Star
## 1         1    *
## ---
## 'CI': 95%-CI for one-sided and 97.5%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 97.5%;
## for two-sided hypotheses, the value tested against lies outside the 97.5%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
# COMP(post-vol) < ADHD(post-vol)
as.data.frame(t(df.des %>% 
    ungroup() %>%
    filter((diagnosis == "ADHD" | diagnosis == "COMP") & 
             ((phase == "volatile") | (phase == "postvolatile"))) %>%
    group_by(diagnosis, phase) %>%
    summarise(across(where(is.numeric), ~ mean(.x))) %>%
    group_by(diagnosis) %>%
    summarise(across(where(is.numeric), ~ diff(.x))) %>% # post - volatile
    select(where(is.numeric)) %>%
    map_df(~ diff(.x)))) %>% # COMP - ADHD
  filter(V1 != 0) 
##                   V1
## diagnosis1:phase1  2
## diagnosis2:phase1  1
## diagnosis1:phase2  4
## diagnosis2:phase2  2
e5 = hypothesis(m.var, "0 > 2*diagnosis1:phase1 + diagnosis2:phase1 +
                + 2*(2*diagnosis1:phase2 + diagnosis2:phase2)", 
               alpha = 0.025)
e5
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (0)-(2*diagnosis1... > 0     0.12      0.07    -0.02     0.25       22.9
##   Post.Prob Star
## 1      0.96     
## ---
## 'CI': 95%-CI for one-sided and 97.5%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 97.5%;
## for two-sided hypotheses, the value tested against lies outside the 97.5%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
# COMP(post-vol) < BOTH(post-vol)
e6 = hypothesis(m.var, "0 > diagnosis1:phase1 + 2*(diagnosis2:phase1 +
                + diagnosis1:phase2 + 2*diagnosis2:phase2)", 
                alpha = 0.025)
e6
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (0)-(diagnosis1:p... > 0     0.05      0.07    -0.08     0.18       3.18
##   Post.Prob Star
## 1      0.76     
## ---
## 'CI': 95%-CI for one-sided and 97.5%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 97.5%;
## for two-sided hypotheses, the value tested against lies outside the 97.5%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
## extract predicted differences in ms instead of log data
df.new = df.var %>% 
  select(diagnosis, phase, expected) %>% 
  distinct() %>%
  mutate(
    condition = paste(diagnosis, phase, expected, sep = "_")
  )
df.ms = as.data.frame(
  fitted(m.var, summary = F, 
               newdata = df.new %>% select(diagnosis, phase, expected), 
               re_formula = NA))
colnames(df.ms) = df.new$condition

# calculate our difference columns
df.ms = df.ms %>%
  mutate(
    COMP          = rowMeans(across(matches("COMP_.*"))),
    ADHD          = rowMeans(across(matches("ADHD_.*"))),
    BOTH          = rowMeans(across(matches("BOTH_.*"))),
    h1a_COMPvADHD = ADHD - COMP,
    e_COMPvBOTH   = BOTH - COMP,
    e_BOTHvADHD   = BOTH - ADHD,
    COMPvol       = rowMeans(across(matches("COMP_volatile.*"))),
    COMPpre       = rowMeans(across(matches("COMP_prevolatile.*"))),
    COMPpost      = rowMeans(across(matches("COMP_postvolatile.*"))),
    ADHDvol       = rowMeans(across(matches("ADHD_volatile.*"))),
    ADHDpre       = rowMeans(across(matches("ADHD_prevolatile.*"))),
    ADHDpost      = rowMeans(across(matches("ADHD_postvolatile.*"))),
    BOTHvol       = rowMeans(across(matches("BOTH_volatile.*"))),
    BOTHpre       = rowMeans(across(matches("BOTH_prevolatile.*"))),
    BOTHpost      = rowMeans(across(matches("BOTH_postvolatile.*"))),
    e_volvpreCOMPvADHD   = (COMPvol  - COMPpre) - (ADHDvol  - ADHDpre),
    e_postvvolCOMPvADHD  = (COMPpost - COMPvol) - (ADHDpost - ADHDvol),
    e_volvpreCOMPvBOTH   = (COMPvol  - COMPpre) - (BOTHvol  - BOTHpre),
    e_postvvolCOMPvBOTH  = (COMPpost - COMPvol) - (BOTHpost - BOTHvol),
    e_volvpreBOTHvADHD   = (BOTHvol  - BOTHpre) - (ADHDvol  - ADHDpre),
    e_postvvolBOTHvADHD  = (BOTHpost - BOTHvol) - (ADHDpost - ADHDvol)
    )

kable(df.ms %>% 
  select(starts_with("e") | starts_with("h")) %>% 
  summarise_all(list(lower.ci = lower_ci, mean = mean, upper.ci = upper_ci)) %>%
  t %>% 
  as.data.frame %>% 
  rownames_to_column() %>%
  separate(rowname, into = c("id", "comparison", "stat"), sep = "_") %>%
  pivot_wider(names_from = stat, values_from = V1) %>%
  arrange(comparison), digits = 2)
id comparison lower.ci mean upper.ci
e BOTHvADHD -32.30 -7.10 18.36
h1a COMPvADHD -0.58 23.43 47.23
e COMPvBOTH -6.40 16.33 40.09
e postvvolBOTHvADHD -22.91 -7.22 8.29
e postvvolCOMPvADHD -26.07 -10.90 3.92
e postvvolCOMPvBOTH -17.74 -3.68 10.42
e volvpreBOTHvADHD -20.00 -4.50 10.95
e volvpreCOMPvADHD -1.24 13.23 27.92
e volvpreCOMPvBOTH 3.57 17.72 31.92
# calculate effect sizes
df.effect = post.draws %>%
  mutate(across(starts_with("sd")|starts_with("sigma"), ~.^2)) %>%
  mutate(
    sumvar = sqrt(rowSums(select(., starts_with("sd")|starts_with("sigma")))),
    h1a    = (2*b_diagnosis1 + b_diagnosis2) / sumvar,
    e1     = (b_diagnosis1 + 2*b_diagnosis2) / sumvar,
    e2     = (b_diagnosis1 - b_diagnosis2) / sumvar,
    e3     = (`b_diagnosis1:phase1` + 2*`b_diagnosis2:phase1` - 
                `b_diagnosis1:phase2` - 2*`b_diagnosis2:phase2`) / sumvar,
    e4     = (2*`b_diagnosis1:phase1` + `b_diagnosis2:phase1` - 
                2*`b_diagnosis1:phase2` - `b_diagnosis2:phase2`) / sumvar
  )

kable(df.effect %>% select(starts_with("e")|starts_with("h")) %>%
        pivot_longer(cols = everything(), values_to = "estimate") %>%
        group_by(name) %>%
        summarise(
          ci.lo = lower_ci(estimate),
          mean  = mean(estimate),
          ci.hi = upper_ci(estimate),
          interpret = interpret_cohens_d(mean)
        ), digits = 3
)
name ci.lo mean ci.hi interpret
e1 -0.137 0.366 0.882 small
e2 -0.356 0.147 0.642 very small
e3 0.099 0.401 0.711 small
e4 0.012 0.313 0.624 small
h1a -0.006 0.513 1.013 medium

h1a: estimate = -0.21 [-0.39, -0.03], posterior probability = 97.41%

e1 COMP < BOTH: estimate = -0.15 [-0.37, 0.06], posterior probability = 92.2%

e3 COMP(vol-pre) > ADHD(vol-pre): estimate = -0.13 [-0.26, -0.01], posterior probability = 97.9%

e4 COMP(vol-pre) > BOTH(vol-pre): estimate = -0.17 [-0.29, -0.04], posterior probability = 99.51%

2.4 Plots

# rain cloud plot
df.var %>%
  mutate(
    group = case_when(
      diagnosis == "BOTH" ~ "ADHD+ASD",
      T ~ diagnosis
    )
  ) %>%
  ggplot(aes(expected, rt.var, fill = group, colour = group)) + #
  geom_rain(rain.side = 'r',
boxplot.args = list(color = "black", outlier.shape = NA, show.legend = FALSE, alpha = .8),
violin.args = list(color = "black", outlier.shape = NA, alpha = .8),
boxplot.args.pos = list(
  position = ggpp::position_dodgenudge(x = 0, width = 0.3), width = 0.3
),
point.args = list(show_guide = FALSE, alpha = .5),
violin.args.pos = list(
  width = 0.6, position = position_nudge(x = 0.16)),
point.args.pos = list(position = ggpp::position_dodgenudge(x = -0.25, width = 0.1))) +
  scale_fill_manual(values = col.grp) +
  scale_color_manual(values = col.grp) +
  #ylim(0, 1) +
  labs(title = "", x = "", y = "variance (ms)") +
  facet_wrap(. ~ phase) +
  theme_bw() + 
  theme(legend.position = "bottom", plot.title = element_blank(), 
        legend.direction = "horizontal", text = element_text(size = 15),
        legend.title = element_blank())

ggsave("plots/FigRTvar.svg", units = "cm", width = 27, height = 9)

# plot transition effects specifically
df.var %>% 
  group_by(diagnosis, phase) %>%
  summarise(
    var.mn = mean(rt.var),
    var.se = sd(rt.var) / sqrt(n())
  ) %>%
  ggplot(aes(y = var.mn, x = phase, 
             group = diagnosis, colour = diagnosis)) +
  geom_line(position = position_dodge(0.4), linewidth = 1) +
  geom_errorbar(aes(ymin = var.mn - var.se, 
                    ymax = var.mn + var.se), linewidth = 1, width = 0.5, 
                position = position_dodge(0.4)) + 
  geom_point(position = position_dodge(0.4), size = 5) + 
  labs (x = "phase", y = "rt var (ms)") + 
  scale_fill_manual(values = col.grp) +
  scale_color_manual(values = col.grp) +
  theme_bw() + 
  theme(legend.position = "bottom", plot.title = element_text(hjust = 0.5))

3 Reaction times in correct trials

In the preregistration, we noted the following population-level effects for the model of the reaction time variances: group, expectancy, phase and difficulty; as well as the group-level predictores subject and trials.

# figure out slopes for subjects
kable(head(df.pal %>% count(subID, expected)))
subID expected n
1 expected 224
1 unexpected 45
10 expected 225
10 unexpected 43
11 expected 225
11 unexpected 46
kable(head(df.pal %>% count(subID, phase)))
subID phase n
1 prevolatile 66
1 volatile 133
1 postvolatile 70
10 prevolatile 67
10 volatile 139
10 postvolatile 62
kable(head(df.pal %>% count(subID, difficulty)))
subID difficulty n
1 easy 95
1 medium 94
1 difficult 80
10 easy 91
10 medium 89
10 difficult 88
kable(head(df.pal %>% count(subID, expected, phase)))
subID expected phase n
1 expected prevolatile 55
1 expected volatile 111
1 expected postvolatile 58
1 unexpected prevolatile 11
1 unexpected volatile 22
1 unexpected postvolatile 12
kable(head(df.pal %>% count(subID, expected, difficulty)))
subID expected difficulty n
1 expected easy 79
1 expected medium 78
1 expected difficult 67
1 unexpected easy 16
1 unexpected medium 16
1 unexpected difficult 13
kable(head(df.pal %>% count(subID, phase, difficulty)))
subID phase difficulty n
1 prevolatile easy 23
1 prevolatile medium 23
1 prevolatile difficult 20
1 volatile easy 48
1 volatile medium 47
1 volatile difficult 38
kable(head(df.pal %>% count(subID, expected, phase, difficulty)))
subID expected phase difficulty n
1 expected prevolatile easy 19
1 expected prevolatile medium 19
1 expected prevolatile difficult 17
1 expected volatile easy 40
1 expected volatile medium 39
1 expected volatile difficult 32
# figure out slopes for trls
kable(head(df.pal %>% count(trl, diagnosis)))
trl diagnosis n
1 ADHD 14
1 BOTH 20
1 COMP 16
2 ADHD 18
2 BOTH 18
2 COMP 21
# set the formula
f.pal = brms::bf(rt.cor ~ diagnosis * expected * phase * difficulty +
                   (expected * phase * difficulty | subID) + (diagnosis | trl))

# set informed priors based on previous results
priors = c(
  # informative priors based Lawson et al. and Schad, Betancourt & Vasishth (2019)
  prior(normal(6.0,   0.3),   class = Intercept),
  prior(normal(0.0,   0.5),   class = sigma),
  prior(normal(0,     0.1),   class = sd),
  prior(lkj(2),               class = cor),
  prior(normal(100, 100.0),   class = ndt), 
  prior(normal(0.00,  0.04),  class = b)
)

3.1 Posterior predictive checks

As the next step, we fit the model, check whether there are divergence or rhat issues, and then check whether the chains have converged.

# fit the final model
m.pal = brm(f.pal, seed = 4466,
            df.pal, prior = priors,
            family = shifted_lognormal,
            iter = iter, warmup = warm,
            backend = "cmdstanr", threads = threading(8),
            file = file.path(brms_dir, "m_pal"),
            save_pars = save_pars(all = TRUE)
            )
rstan::check_hmc_diagnostics(m.pal$fit)
## 
## Divergences:
## 0 of 18000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 18000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.
# check that rhats are below 1.01
sum(brms::rhat(m.pal) >= 1.01, na.rm = T)
## [1] 0
# check the trace plots
post.draws = as_draws_df(m.pal)
mcmc_trace(post.draws, regex_pars = "^b_",
           facet_args = list(ncol = 6)) +
  scale_x_continuous(breaks=scales::pretty_breaks(n = 3)) +
  scale_y_continuous(breaks=scales::pretty_breaks(n = 3))
## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.

This model has no pathological behaviour with E-BFMI, no divergent sample and no rhats that are higher or equal to 1.01. Therefore, we go ahead and perform our posterior predictive checks.

# get posterior predictions
post.pred = posterior_predict(m.pal, ndraws = nsim)

# check the fit of the predicted data compared to the real data
p1 = pp_check(m.pal, ndraws = nsim) + 
  theme_bw() + theme(legend.position = "none") + xlim(0, 1500)

# distributions of means compared to the real values per group or conditions
p2 = ppc_stat_grouped(df.pal$rt.cor, 
                      post.pred, 
                      df.pal$diagnosis) + 
  theme_bw() + theme(legend.position = "none")
p3 = ppc_stat_grouped(df.pal$rt.cor, 
                      post.pred, 
                      df.pal$expected) + 
  theme_bw() + theme(legend.position = "none")
p4 = ppc_stat_grouped(df.pal$rt.cor, 
                      post.pred, 
                      df.pal$phase) + 
  theme_bw() + theme(legend.position = "none")
p5 = ppc_stat_grouped(df.pal$rt.cor, 
                      post.pred, 
                      df.pal$difficulty) + 
  theme_bw() + theme(legend.position = "none")

p = ggarrange(p1, p2, p3, p4, p5, 
          nrow = 5, ncol = 1, labels = "AUTO")
annotate_figure(p, top = text_grob("Posterior predictive checks: reaction times", 
                                   face = "bold", size = 14))

The simulated data based on the model fits well with the real data, although there are slight deviations specifically for the predictor expectedness. However, we judge this to be acceptable.

3.2 Inferences

Now that we are convinced that we can trust our model, we have a look at its estimate and use the hypothesis function to assess our hypotheses and perform explorative tests.

# print a summary
summary(m.pal)
##  Family: shifted_lognormal 
##   Links: mu = identity; sigma = identity; ndt = identity 
## Formula: rt.cor ~ diagnosis * expected * phase * difficulty + (expected * phase * difficulty | subID) + (diagnosis | trl) 
##    Data: df.pal (Number of observations: 16886) 
##   Draws: 4 chains, each with iter = 6000; warmup = 1500; thin = 1;
##          total post-warmup draws = 18000
## 
## Multilevel Hyperparameters:
## ~subID (Number of levels: 66) 
##                                                                Estimate
## sd(Intercept)                                                      0.19
## sd(expected1)                                                      0.02
## sd(phase1)                                                         0.06
## sd(phase2)                                                         0.02
## sd(difficulty1)                                                    0.01
## sd(difficulty2)                                                    0.01
## sd(expected1:phase1)                                               0.02
## sd(expected1:phase2)                                               0.01
## sd(expected1:difficulty1)                                          0.00
## sd(expected1:difficulty2)                                          0.01
## sd(phase1:difficulty1)                                             0.01
## sd(phase2:difficulty1)                                             0.01
## sd(phase1:difficulty2)                                             0.01
## sd(phase2:difficulty2)                                             0.01
## sd(expected1:phase1:difficulty1)                                   0.01
## sd(expected1:phase2:difficulty1)                                   0.01
## sd(expected1:phase1:difficulty2)                                   0.01
## sd(expected1:phase2:difficulty2)                                   0.01
## cor(Intercept,expected1)                                           0.23
## cor(Intercept,phase1)                                              0.08
## cor(expected1,phase1)                                              0.22
## cor(Intercept,phase2)                                              0.09
## cor(expected1,phase2)                                              0.04
## cor(phase1,phase2)                                                -0.25
## cor(Intercept,difficulty1)                                        -0.25
## cor(expected1,difficulty1)                                        -0.08
## cor(phase1,difficulty1)                                           -0.17
## cor(phase2,difficulty1)                                           -0.07
## cor(Intercept,difficulty2)                                        -0.04
## cor(expected1,difficulty2)                                        -0.11
## cor(phase1,difficulty2)                                           -0.17
## cor(phase2,difficulty2)                                            0.09
## cor(difficulty1,difficulty2)                                       0.02
## cor(Intercept,expected1:phase1)                                    0.03
## cor(expected1,expected1:phase1)                                    0.17
## cor(phase1,expected1:phase1)                                       0.09
## cor(phase2,expected1:phase1)                                      -0.17
## cor(difficulty1,expected1:phase1)                                  0.06
## cor(difficulty2,expected1:phase1)                                 -0.02
## cor(Intercept,expected1:phase2)                                    0.21
## cor(expected1,expected1:phase2)                                    0.07
## cor(phase1,expected1:phase2)                                      -0.05
## cor(phase2,expected1:phase2)                                       0.04
## cor(difficulty1,expected1:phase2)                                  0.05
## cor(difficulty2,expected1:phase2)                                  0.13
## cor(expected1:phase1,expected1:phase2)                             0.01
## cor(Intercept,expected1:difficulty1)                              -0.07
## cor(expected1,expected1:difficulty1)                               0.03
## cor(phase1,expected1:difficulty1)                                  0.03
## cor(phase2,expected1:difficulty1)                                  0.00
## cor(difficulty1,expected1:difficulty1)                            -0.01
## cor(difficulty2,expected1:difficulty1)                            -0.00
## cor(expected1:phase1,expected1:difficulty1)                        0.04
## cor(expected1:phase2,expected1:difficulty1)                        0.01
## cor(Intercept,expected1:difficulty2)                               0.01
## cor(expected1,expected1:difficulty2)                               0.05
## cor(phase1,expected1:difficulty2)                                  0.10
## cor(phase2,expected1:difficulty2)                                  0.07
## cor(difficulty1,expected1:difficulty2)                             0.02
## cor(difficulty2,expected1:difficulty2)                            -0.01
## cor(expected1:phase1,expected1:difficulty2)                        0.06
## cor(expected1:phase2,expected1:difficulty2)                        0.07
## cor(expected1:difficulty1,expected1:difficulty2)                  -0.01
## cor(Intercept,phase1:difficulty1)                                 -0.09
## cor(expected1,phase1:difficulty1)                                  0.06
## cor(phase1,phase1:difficulty1)                                     0.15
## cor(phase2,phase1:difficulty1)                                    -0.11
## cor(difficulty1,phase1:difficulty1)                                0.07
## cor(difficulty2,phase1:difficulty1)                               -0.06
## cor(expected1:phase1,phase1:difficulty1)                           0.07
## cor(expected1:phase2,phase1:difficulty1)                          -0.04
## cor(expected1:difficulty1,phase1:difficulty1)                      0.05
## cor(expected1:difficulty2,phase1:difficulty1)                      0.01
## cor(Intercept,phase2:difficulty1)                                 -0.05
## cor(expected1,phase2:difficulty1)                                  0.05
## cor(phase1,phase2:difficulty1)                                    -0.11
## cor(phase2,phase2:difficulty1)                                     0.03
## cor(difficulty1,phase2:difficulty1)                                0.09
## cor(difficulty2,phase2:difficulty1)                                0.01
## cor(expected1:phase1,phase2:difficulty1)                           0.01
## cor(expected1:phase2,phase2:difficulty1)                           0.08
## cor(expected1:difficulty1,phase2:difficulty1)                      0.03
## cor(expected1:difficulty2,phase2:difficulty1)                     -0.03
## cor(phase1:difficulty1,phase2:difficulty1)                        -0.02
## cor(Intercept,phase1:difficulty2)                                  0.06
## cor(expected1,phase1:difficulty2)                                 -0.03
## cor(phase1,phase1:difficulty2)                                     0.04
## cor(phase2,phase1:difficulty2)                                    -0.06
## cor(difficulty1,phase1:difficulty2)                               -0.05
## cor(difficulty2,phase1:difficulty2)                                0.01
## cor(expected1:phase1,phase1:difficulty2)                          -0.03
## cor(expected1:phase2,phase1:difficulty2)                          -0.03
## cor(expected1:difficulty1,phase1:difficulty2)                     -0.04
## cor(expected1:difficulty2,phase1:difficulty2)                     -0.00
## cor(phase1:difficulty1,phase1:difficulty2)                        -0.05
## cor(phase2:difficulty1,phase1:difficulty2)                        -0.05
## cor(Intercept,phase2:difficulty2)                                  0.23
## cor(expected1,phase2:difficulty2)                                  0.13
## cor(phase1,phase2:difficulty2)                                    -0.08
## cor(phase2,phase2:difficulty2)                                     0.04
## cor(difficulty1,phase2:difficulty2)                                0.04
## cor(difficulty2,phase2:difficulty2)                                0.02
## cor(expected1:phase1,phase2:difficulty2)                           0.09
## cor(expected1:phase2,phase2:difficulty2)                           0.19
## cor(expected1:difficulty1,phase2:difficulty2)                     -0.02
## cor(expected1:difficulty2,phase2:difficulty2)                      0.01
## cor(phase1:difficulty1,phase2:difficulty2)                        -0.09
## cor(phase2:difficulty1,phase2:difficulty2)                        -0.04
## cor(phase1:difficulty2,phase2:difficulty2)                        -0.05
## cor(Intercept,expected1:phase1:difficulty1)                        0.06
## cor(expected1,expected1:phase1:difficulty1)                        0.07
## cor(phase1,expected1:phase1:difficulty1)                           0.13
## cor(phase2,expected1:phase1:difficulty1)                          -0.01
## cor(difficulty1,expected1:phase1:difficulty1)                     -0.01
## cor(difficulty2,expected1:phase1:difficulty1)                     -0.01
## cor(expected1:phase1,expected1:phase1:difficulty1)                 0.03
## cor(expected1:phase2,expected1:phase1:difficulty1)                 0.02
## cor(expected1:difficulty1,expected1:phase1:difficulty1)            0.02
## cor(expected1:difficulty2,expected1:phase1:difficulty1)            0.04
## cor(phase1:difficulty1,expected1:phase1:difficulty1)              -0.05
## cor(phase2:difficulty1,expected1:phase1:difficulty1)              -0.04
## cor(phase1:difficulty2,expected1:phase1:difficulty1)              -0.02
## cor(phase2:difficulty2,expected1:phase1:difficulty1)              -0.00
## cor(Intercept,expected1:phase2:difficulty1)                        0.04
## cor(expected1,expected1:phase2:difficulty1)                        0.10
## cor(phase1,expected1:phase2:difficulty1)                           0.06
## cor(phase2,expected1:phase2:difficulty1)                           0.05
## cor(difficulty1,expected1:phase2:difficulty1)                     -0.04
## cor(difficulty2,expected1:phase2:difficulty1)                     -0.03
## cor(expected1:phase1,expected1:phase2:difficulty1)                -0.01
## cor(expected1:phase2,expected1:phase2:difficulty1)                -0.01
## cor(expected1:difficulty1,expected1:phase2:difficulty1)           -0.02
## cor(expected1:difficulty2,expected1:phase2:difficulty1)            0.01
## cor(phase1:difficulty1,expected1:phase2:difficulty1)              -0.02
## cor(phase2:difficulty1,expected1:phase2:difficulty1)              -0.05
## cor(phase1:difficulty2,expected1:phase2:difficulty1)              -0.01
## cor(phase2:difficulty2,expected1:phase2:difficulty1)              -0.01
## cor(expected1:phase1:difficulty1,expected1:phase2:difficulty1)    -0.00
## cor(Intercept,expected1:phase1:difficulty2)                        0.02
## cor(expected1,expected1:phase1:difficulty2)                       -0.04
## cor(phase1,expected1:phase1:difficulty2)                           0.02
## cor(phase2,expected1:phase1:difficulty2)                          -0.01
## cor(difficulty1,expected1:phase1:difficulty2)                     -0.02
## cor(difficulty2,expected1:phase1:difficulty2)                      0.03
## cor(expected1:phase1,expected1:phase1:difficulty2)                -0.03
## cor(expected1:phase2,expected1:phase1:difficulty2)                -0.03
## cor(expected1:difficulty1,expected1:phase1:difficulty2)           -0.02
## cor(expected1:difficulty2,expected1:phase1:difficulty2)            0.01
## cor(phase1:difficulty1,expected1:phase1:difficulty2)              -0.04
## cor(phase2:difficulty1,expected1:phase1:difficulty2)              -0.03
## cor(phase1:difficulty2,expected1:phase1:difficulty2)              -0.01
## cor(phase2:difficulty2,expected1:phase1:difficulty2)              -0.03
## cor(expected1:phase1:difficulty1,expected1:phase1:difficulty2)    -0.03
## cor(expected1:phase2:difficulty1,expected1:phase1:difficulty2)    -0.01
## cor(Intercept,expected1:phase2:difficulty2)                        0.04
## cor(expected1,expected1:phase2:difficulty2)                        0.07
## cor(phase1,expected1:phase2:difficulty2)                           0.09
## cor(phase2,expected1:phase2:difficulty2)                           0.06
## cor(difficulty1,expected1:phase2:difficulty2)                     -0.03
## cor(difficulty2,expected1:phase2:difficulty2)                     -0.04
## cor(expected1:phase1,expected1:phase2:difficulty2)                 0.01
## cor(expected1:phase2,expected1:phase2:difficulty2)                 0.01
## cor(expected1:difficulty1,expected1:phase2:difficulty2)           -0.01
## cor(expected1:difficulty2,expected1:phase2:difficulty2)           -0.00
## cor(phase1:difficulty1,expected1:phase2:difficulty2)              -0.00
## cor(phase2:difficulty1,expected1:phase2:difficulty2)              -0.04
## cor(phase1:difficulty2,expected1:phase2:difficulty2)              -0.03
## cor(phase2:difficulty2,expected1:phase2:difficulty2)              -0.07
## cor(expected1:phase1:difficulty1,expected1:phase2:difficulty2)     0.01
## cor(expected1:phase2:difficulty1,expected1:phase2:difficulty2)    -0.01
## cor(expected1:phase1:difficulty2,expected1:phase2:difficulty2)    -0.03
##                                                                Est.Error
## sd(Intercept)                                                       0.02
## sd(expected1)                                                       0.00
## sd(phase1)                                                          0.01
## sd(phase2)                                                          0.00
## sd(difficulty1)                                                     0.00
## sd(difficulty2)                                                     0.00
## sd(expected1:phase1)                                                0.01
## sd(expected1:phase2)                                                0.01
## sd(expected1:difficulty1)                                           0.00
## sd(expected1:difficulty2)                                           0.00
## sd(phase1:difficulty1)                                              0.01
## sd(phase2:difficulty1)                                              0.00
## sd(phase1:difficulty2)                                              0.00
## sd(phase2:difficulty2)                                              0.01
## sd(expected1:phase1:difficulty1)                                    0.00
## sd(expected1:phase2:difficulty1)                                    0.00
## sd(expected1:phase1:difficulty2)                                    0.00
## sd(expected1:phase2:difficulty2)                                    0.00
## cor(Intercept,expected1)                                            0.14
## cor(Intercept,phase1)                                               0.12
## cor(expected1,phase1)                                               0.15
## cor(Intercept,phase2)                                               0.14
## cor(expected1,phase2)                                               0.18
## cor(phase1,phase2)                                                  0.15
## cor(Intercept,difficulty1)                                          0.17
## cor(expected1,difficulty1)                                          0.19
## cor(phase1,difficulty1)                                             0.17
## cor(phase2,difficulty1)                                             0.19
## cor(Intercept,difficulty2)                                          0.19
## cor(expected1,difficulty2)                                          0.20
## cor(phase1,difficulty2)                                             0.19
## cor(phase2,difficulty2)                                             0.20
## cor(difficulty1,difficulty2)                                        0.21
## cor(Intercept,expected1:phase1)                                     0.17
## cor(expected1,expected1:phase1)                                     0.19
## cor(phase1,expected1:phase1)                                        0.18
## cor(phase2,expected1:phase1)                                        0.19
## cor(difficulty1,expected1:phase1)                                   0.20
## cor(difficulty2,expected1:phase1)                                   0.21
## cor(Intercept,expected1:phase2)                                     0.17
## cor(expected1,expected1:phase2)                                     0.19
## cor(phase1,expected1:phase2)                                        0.17
## cor(phase2,expected1:phase2)                                        0.19
## cor(difficulty1,expected1:phase2)                                   0.20
## cor(difficulty2,expected1:phase2)                                   0.21
## cor(expected1:phase1,expected1:phase2)                              0.20
## cor(Intercept,expected1:difficulty1)                                0.21
## cor(expected1,expected1:difficulty1)                                0.22
## cor(phase1,expected1:difficulty1)                                   0.21
## cor(phase2,expected1:difficulty1)                                   0.21
## cor(difficulty1,expected1:difficulty1)                              0.22
## cor(difficulty2,expected1:difficulty1)                              0.21
## cor(expected1:phase1,expected1:difficulty1)                         0.22
## cor(expected1:phase2,expected1:difficulty1)                         0.21
## cor(Intercept,expected1:difficulty2)                                0.19
## cor(expected1,expected1:difficulty2)                                0.21
## cor(phase1,expected1:difficulty2)                                   0.20
## cor(phase2,expected1:difficulty2)                                   0.21
## cor(difficulty1,expected1:difficulty2)                              0.21
## cor(difficulty2,expected1:difficulty2)                              0.22
## cor(expected1:phase1,expected1:difficulty2)                         0.21
## cor(expected1:phase2,expected1:difficulty2)                         0.21
## cor(expected1:difficulty1,expected1:difficulty2)                    0.22
## cor(Intercept,phase1:difficulty1)                                   0.20
## cor(expected1,phase1:difficulty1)                                   0.21
## cor(phase1,phase1:difficulty1)                                      0.21
## cor(phase2,phase1:difficulty1)                                      0.21
## cor(difficulty1,phase1:difficulty1)                                 0.21
## cor(difficulty2,phase1:difficulty1)                                 0.21
## cor(expected1:phase1,phase1:difficulty1)                            0.21
## cor(expected1:phase2,phase1:difficulty1)                            0.21
## cor(expected1:difficulty1,phase1:difficulty1)                       0.22
## cor(expected1:difficulty2,phase1:difficulty1)                       0.22
## cor(Intercept,phase2:difficulty1)                                   0.20
## cor(expected1,phase2:difficulty1)                                   0.21
## cor(phase1,phase2:difficulty1)                                      0.20
## cor(phase2,phase2:difficulty1)                                      0.21
## cor(difficulty1,phase2:difficulty1)                                 0.21
## cor(difficulty2,phase2:difficulty1)                                 0.21
## cor(expected1:phase1,phase2:difficulty1)                            0.21
## cor(expected1:phase2,phase2:difficulty1)                            0.21
## cor(expected1:difficulty1,phase2:difficulty1)                       0.22
## cor(expected1:difficulty2,phase2:difficulty1)                       0.22
## cor(phase1:difficulty1,phase2:difficulty1)                          0.22
## cor(Intercept,phase1:difficulty2)                                   0.21
## cor(expected1,phase1:difficulty2)                                   0.21
## cor(phase1,phase1:difficulty2)                                      0.21
## cor(phase2,phase1:difficulty2)                                      0.21
## cor(difficulty1,phase1:difficulty2)                                 0.22
## cor(difficulty2,phase1:difficulty2)                                 0.22
## cor(expected1:phase1,phase1:difficulty2)                            0.22
## cor(expected1:phase2,phase1:difficulty2)                            0.22
## cor(expected1:difficulty1,phase1:difficulty2)                       0.22
## cor(expected1:difficulty2,phase1:difficulty2)                       0.22
## cor(phase1:difficulty1,phase1:difficulty2)                          0.22
## cor(phase2:difficulty1,phase1:difficulty2)                          0.22
## cor(Intercept,phase2:difficulty2)                                   0.17
## cor(expected1,phase2:difficulty2)                                   0.19
## cor(phase1,phase2:difficulty2)                                      0.18
## cor(phase2,phase2:difficulty2)                                      0.19
## cor(difficulty1,phase2:difficulty2)                                 0.20
## cor(difficulty2,phase2:difficulty2)                                 0.21
## cor(expected1:phase1,phase2:difficulty2)                            0.20
## cor(expected1:phase2,phase2:difficulty2)                            0.20
## cor(expected1:difficulty1,phase2:difficulty2)                       0.22
## cor(expected1:difficulty2,phase2:difficulty2)                       0.21
## cor(phase1:difficulty1,phase2:difficulty2)                          0.22
## cor(phase2:difficulty1,phase2:difficulty2)                          0.21
## cor(phase1:difficulty2,phase2:difficulty2)                          0.22
## cor(Intercept,expected1:phase1:difficulty1)                         0.21
## cor(expected1,expected1:phase1:difficulty1)                         0.22
## cor(phase1,expected1:phase1:difficulty1)                            0.22
## cor(phase2,expected1:phase1:difficulty1)                            0.22
## cor(difficulty1,expected1:phase1:difficulty1)                       0.21
## cor(difficulty2,expected1:phase1:difficulty1)                       0.22
## cor(expected1:phase1,expected1:phase1:difficulty1)                  0.22
## cor(expected1:phase2,expected1:phase1:difficulty1)                  0.21
## cor(expected1:difficulty1,expected1:phase1:difficulty1)             0.22
## cor(expected1:difficulty2,expected1:phase1:difficulty1)             0.22
## cor(phase1:difficulty1,expected1:phase1:difficulty1)                0.22
## cor(phase2:difficulty1,expected1:phase1:difficulty1)                0.22
## cor(phase1:difficulty2,expected1:phase1:difficulty1)                0.22
## cor(phase2:difficulty2,expected1:phase1:difficulty1)                0.22
## cor(Intercept,expected1:phase2:difficulty1)                         0.21
## cor(expected1,expected1:phase2:difficulty1)                         0.22
## cor(phase1,expected1:phase2:difficulty1)                            0.21
## cor(phase2,expected1:phase2:difficulty1)                            0.21
## cor(difficulty1,expected1:phase2:difficulty1)                       0.22
## cor(difficulty2,expected1:phase2:difficulty1)                       0.22
## cor(expected1:phase1,expected1:phase2:difficulty1)                  0.22
## cor(expected1:phase2,expected1:phase2:difficulty1)                  0.22
## cor(expected1:difficulty1,expected1:phase2:difficulty1)             0.22
## cor(expected1:difficulty2,expected1:phase2:difficulty1)             0.22
## cor(phase1:difficulty1,expected1:phase2:difficulty1)                0.22
## cor(phase2:difficulty1,expected1:phase2:difficulty1)                0.22
## cor(phase1:difficulty2,expected1:phase2:difficulty1)                0.22
## cor(phase2:difficulty2,expected1:phase2:difficulty1)                0.22
## cor(expected1:phase1:difficulty1,expected1:phase2:difficulty1)      0.22
## cor(Intercept,expected1:phase1:difficulty2)                         0.21
## cor(expected1,expected1:phase1:difficulty2)                         0.21
## cor(phase1,expected1:phase1:difficulty2)                            0.21
## cor(phase2,expected1:phase1:difficulty2)                            0.21
## cor(difficulty1,expected1:phase1:difficulty2)                       0.21
## cor(difficulty2,expected1:phase1:difficulty2)                       0.22
## cor(expected1:phase1,expected1:phase1:difficulty2)                  0.21
## cor(expected1:phase2,expected1:phase1:difficulty2)                  0.22
## cor(expected1:difficulty1,expected1:phase1:difficulty2)             0.22
## cor(expected1:difficulty2,expected1:phase1:difficulty2)             0.22
## cor(phase1:difficulty1,expected1:phase1:difficulty2)                0.22
## cor(phase2:difficulty1,expected1:phase1:difficulty2)                0.22
## cor(phase1:difficulty2,expected1:phase1:difficulty2)                0.22
## cor(phase2:difficulty2,expected1:phase1:difficulty2)                0.22
## cor(expected1:phase1:difficulty1,expected1:phase1:difficulty2)      0.22
## cor(expected1:phase2:difficulty1,expected1:phase1:difficulty2)      0.21
## cor(Intercept,expected1:phase2:difficulty2)                         0.21
## cor(expected1,expected1:phase2:difficulty2)                         0.21
## cor(phase1,expected1:phase2:difficulty2)                            0.21
## cor(phase2,expected1:phase2:difficulty2)                            0.21
## cor(difficulty1,expected1:phase2:difficulty2)                       0.21
## cor(difficulty2,expected1:phase2:difficulty2)                       0.22
## cor(expected1:phase1,expected1:phase2:difficulty2)                  0.21
## cor(expected1:phase2,expected1:phase2:difficulty2)                  0.21
## cor(expected1:difficulty1,expected1:phase2:difficulty2)             0.22
## cor(expected1:difficulty2,expected1:phase2:difficulty2)             0.22
## cor(phase1:difficulty1,expected1:phase2:difficulty2)                0.22
## cor(phase2:difficulty1,expected1:phase2:difficulty2)                0.22
## cor(phase1:difficulty2,expected1:phase2:difficulty2)                0.22
## cor(phase2:difficulty2,expected1:phase2:difficulty2)                0.22
## cor(expected1:phase1:difficulty1,expected1:phase2:difficulty2)      0.22
## cor(expected1:phase2:difficulty1,expected1:phase2:difficulty2)      0.22
## cor(expected1:phase1:difficulty2,expected1:phase2:difficulty2)      0.22
##                                                                l-95% CI
## sd(Intercept)                                                      0.16
## sd(expected1)                                                      0.01
## sd(phase1)                                                         0.05
## sd(phase2)                                                         0.01
## sd(difficulty1)                                                    0.00
## sd(difficulty2)                                                    0.00
## sd(expected1:phase1)                                               0.00
## sd(expected1:phase2)                                               0.00
## sd(expected1:difficulty1)                                          0.00
## sd(expected1:difficulty2)                                          0.00
## sd(phase1:difficulty1)                                             0.00
## sd(phase2:difficulty1)                                             0.00
## sd(phase1:difficulty2)                                             0.00
## sd(phase2:difficulty2)                                             0.00
## sd(expected1:phase1:difficulty1)                                   0.00
## sd(expected1:phase2:difficulty1)                                   0.00
## sd(expected1:phase1:difficulty2)                                   0.00
## sd(expected1:phase2:difficulty2)                                   0.00
## cor(Intercept,expected1)                                          -0.06
## cor(Intercept,phase1)                                             -0.15
## cor(expected1,phase1)                                             -0.09
## cor(Intercept,phase2)                                             -0.19
## cor(expected1,phase2)                                             -0.31
## cor(phase1,phase2)                                                -0.52
## cor(Intercept,difficulty1)                                        -0.55
## cor(expected1,difficulty1)                                        -0.44
## cor(phase1,difficulty1)                                           -0.50
## cor(phase2,difficulty1)                                           -0.44
## cor(Intercept,difficulty2)                                        -0.40
## cor(expected1,difficulty2)                                        -0.49
## cor(phase1,difficulty2)                                           -0.52
## cor(phase2,difficulty2)                                           -0.31
## cor(difficulty1,difficulty2)                                      -0.38
## cor(Intercept,expected1:phase1)                                   -0.30
## cor(expected1,expected1:phase1)                                   -0.22
## cor(phase1,expected1:phase1)                                      -0.25
## cor(phase2,expected1:phase1)                                      -0.53
## cor(difficulty1,expected1:phase1)                                 -0.33
## cor(difficulty2,expected1:phase1)                                 -0.42
## cor(Intercept,expected1:phase2)                                   -0.14
## cor(expected1,expected1:phase2)                                   -0.30
## cor(phase1,expected1:phase2)                                      -0.38
## cor(phase2,expected1:phase2)                                      -0.33
## cor(difficulty1,expected1:phase2)                                 -0.35
## cor(difficulty2,expected1:phase2)                                 -0.29
## cor(expected1:phase1,expected1:phase2)                            -0.38
## cor(Intercept,expected1:difficulty1)                              -0.47
## cor(expected1,expected1:difficulty1)                              -0.39
## cor(phase1,expected1:difficulty1)                                 -0.40
## cor(phase2,expected1:difficulty1)                                 -0.42
## cor(difficulty1,expected1:difficulty1)                            -0.43
## cor(difficulty2,expected1:difficulty1)                            -0.42
## cor(expected1:phase1,expected1:difficulty1)                       -0.39
## cor(expected1:phase2,expected1:difficulty1)                       -0.40
## cor(Intercept,expected1:difficulty2)                              -0.37
## cor(expected1,expected1:difficulty2)                              -0.36
## cor(phase1,expected1:difficulty2)                                 -0.31
## cor(phase2,expected1:difficulty2)                                 -0.34
## cor(difficulty1,expected1:difficulty2)                            -0.39
## cor(difficulty2,expected1:difficulty2)                            -0.42
## cor(expected1:phase1,expected1:difficulty2)                       -0.36
## cor(expected1:phase2,expected1:difficulty2)                       -0.34
## cor(expected1:difficulty1,expected1:difficulty2)                  -0.42
## cor(Intercept,phase1:difficulty1)                                 -0.46
## cor(expected1,phase1:difficulty1)                                 -0.36
## cor(phase1,phase1:difficulty1)                                    -0.28
## cor(phase2,phase1:difficulty1)                                    -0.51
## cor(difficulty1,phase1:difficulty1)                               -0.35
## cor(difficulty2,phase1:difficulty1)                               -0.48
## cor(expected1:phase1,phase1:difficulty1)                          -0.35
## cor(expected1:phase2,phase1:difficulty1)                          -0.44
## cor(expected1:difficulty1,phase1:difficulty1)                     -0.38
## cor(expected1:difficulty2,phase1:difficulty1)                     -0.41
## cor(Intercept,phase2:difficulty1)                                 -0.42
## cor(expected1,phase2:difficulty1)                                 -0.37
## cor(phase1,phase2:difficulty1)                                    -0.49
## cor(phase2,phase2:difficulty1)                                    -0.37
## cor(difficulty1,phase2:difficulty1)                               -0.34
## cor(difficulty2,phase2:difficulty1)                               -0.40
## cor(expected1:phase1,phase2:difficulty1)                          -0.40
## cor(expected1:phase2,phase2:difficulty1)                          -0.35
## cor(expected1:difficulty1,phase2:difficulty1)                     -0.40
## cor(expected1:difficulty2,phase2:difficulty1)                     -0.45
## cor(phase1:difficulty1,phase2:difficulty1)                        -0.44
## cor(Intercept,phase1:difficulty2)                                 -0.36
## cor(expected1,phase1:difficulty2)                                 -0.43
## cor(phase1,phase1:difficulty2)                                    -0.37
## cor(phase2,phase1:difficulty2)                                    -0.46
## cor(difficulty1,phase1:difficulty2)                               -0.46
## cor(difficulty2,phase1:difficulty2)                               -0.41
## cor(expected1:phase1,phase1:difficulty2)                          -0.45
## cor(expected1:phase2,phase1:difficulty2)                          -0.44
## cor(expected1:difficulty1,phase1:difficulty2)                     -0.46
## cor(expected1:difficulty2,phase1:difficulty2)                     -0.42
## cor(phase1:difficulty1,phase1:difficulty2)                        -0.48
## cor(phase2:difficulty1,phase1:difficulty2)                        -0.47
## cor(Intercept,phase2:difficulty2)                                 -0.14
## cor(expected1,phase2:difficulty2)                                 -0.26
## cor(phase1,phase2:difficulty2)                                    -0.42
## cor(phase2,phase2:difficulty2)                                    -0.33
## cor(difficulty1,phase2:difficulty2)                               -0.36
## cor(difficulty2,phase2:difficulty2)                               -0.38
## cor(expected1:phase1,phase2:difficulty2)                          -0.32
## cor(expected1:phase2,phase2:difficulty2)                          -0.22
## cor(expected1:difficulty1,phase2:difficulty2)                     -0.43
## cor(expected1:difficulty2,phase2:difficulty2)                     -0.41
## cor(phase1:difficulty1,phase2:difficulty2)                        -0.49
## cor(phase2:difficulty1,phase2:difficulty2)                        -0.44
## cor(phase1:difficulty2,phase2:difficulty2)                        -0.46
## cor(Intercept,expected1:phase1:difficulty1)                       -0.36
## cor(expected1,expected1:phase1:difficulty1)                       -0.37
## cor(phase1,expected1:phase1:difficulty1)                          -0.33
## cor(phase2,expected1:phase1:difficulty1)                          -0.43
## cor(difficulty1,expected1:phase1:difficulty1)                     -0.42
## cor(difficulty2,expected1:phase1:difficulty1)                     -0.43
## cor(expected1:phase1,expected1:phase1:difficulty1)                -0.40
## cor(expected1:phase2,expected1:phase1:difficulty1)                -0.39
## cor(expected1:difficulty1,expected1:phase1:difficulty1)           -0.41
## cor(expected1:difficulty2,expected1:phase1:difficulty1)           -0.40
## cor(phase1:difficulty1,expected1:phase1:difficulty1)              -0.47
## cor(phase2:difficulty1,expected1:phase1:difficulty1)              -0.46
## cor(phase1:difficulty2,expected1:phase1:difficulty1)              -0.44
## cor(phase2:difficulty2,expected1:phase1:difficulty1)              -0.41
## cor(Intercept,expected1:phase2:difficulty1)                       -0.37
## cor(expected1,expected1:phase2:difficulty1)                       -0.34
## cor(phase1,expected1:phase2:difficulty1)                          -0.37
## cor(phase2,expected1:phase2:difficulty1)                          -0.37
## cor(difficulty1,expected1:phase2:difficulty1)                     -0.45
## cor(difficulty2,expected1:phase2:difficulty1)                     -0.44
## cor(expected1:phase1,expected1:phase2:difficulty1)                -0.44
## cor(expected1:phase2,expected1:phase2:difficulty1)                -0.42
## cor(expected1:difficulty1,expected1:phase2:difficulty1)           -0.43
## cor(expected1:difficulty2,expected1:phase2:difficulty1)           -0.41
## cor(phase1:difficulty1,expected1:phase2:difficulty1)              -0.44
## cor(phase2:difficulty1,expected1:phase2:difficulty1)              -0.47
## cor(phase1:difficulty2,expected1:phase2:difficulty1)              -0.44
## cor(phase2:difficulty2,expected1:phase2:difficulty1)              -0.42
## cor(expected1:phase1:difficulty1,expected1:phase2:difficulty1)    -0.42
## cor(Intercept,expected1:phase1:difficulty2)                       -0.38
## cor(expected1,expected1:phase1:difficulty2)                       -0.45
## cor(phase1,expected1:phase1:difficulty2)                          -0.39
## cor(phase2,expected1:phase1:difficulty2)                          -0.43
## cor(difficulty1,expected1:phase1:difficulty2)                     -0.43
## cor(difficulty2,expected1:phase1:difficulty2)                     -0.40
## cor(expected1:phase1,expected1:phase1:difficulty2)                -0.44
## cor(expected1:phase2,expected1:phase1:difficulty2)                -0.45
## cor(expected1:difficulty1,expected1:phase1:difficulty2)           -0.43
## cor(expected1:difficulty2,expected1:phase1:difficulty2)           -0.41
## cor(phase1:difficulty1,expected1:phase1:difficulty2)              -0.46
## cor(phase2:difficulty1,expected1:phase1:difficulty2)              -0.45
## cor(phase1:difficulty2,expected1:phase1:difficulty2)              -0.44
## cor(phase2:difficulty2,expected1:phase1:difficulty2)              -0.45
## cor(expected1:phase1:difficulty1,expected1:phase1:difficulty2)    -0.46
## cor(expected1:phase2:difficulty1,expected1:phase1:difficulty2)    -0.43
## cor(Intercept,expected1:phase2:difficulty2)                       -0.37
## cor(expected1,expected1:phase2:difficulty2)                       -0.35
## cor(phase1,expected1:phase2:difficulty2)                          -0.33
## cor(phase2,expected1:phase2:difficulty2)                          -0.36
## cor(difficulty1,expected1:phase2:difficulty2)                     -0.44
## cor(difficulty2,expected1:phase2:difficulty2)                     -0.46
## cor(expected1:phase1,expected1:phase2:difficulty2)                -0.40
## cor(expected1:phase2,expected1:phase2:difficulty2)                -0.42
## cor(expected1:difficulty1,expected1:phase2:difficulty2)           -0.42
## cor(expected1:difficulty2,expected1:phase2:difficulty2)           -0.42
## cor(phase1:difficulty1,expected1:phase2:difficulty2)              -0.42
## cor(phase2:difficulty1,expected1:phase2:difficulty2)              -0.46
## cor(phase1:difficulty2,expected1:phase2:difficulty2)              -0.45
## cor(phase2:difficulty2,expected1:phase2:difficulty2)              -0.49
## cor(expected1:phase1:difficulty1,expected1:phase2:difficulty2)    -0.41
## cor(expected1:phase2:difficulty1,expected1:phase2:difficulty2)    -0.42
## cor(expected1:phase1:difficulty2,expected1:phase2:difficulty2)    -0.46
##                                                                u-95% CI Rhat
## sd(Intercept)                                                      0.23 1.00
## sd(expected1)                                                      0.03 1.00
## sd(phase1)                                                         0.07 1.00
## sd(phase2)                                                         0.03 1.00
## sd(difficulty1)                                                    0.02 1.00
## sd(difficulty2)                                                    0.02 1.00
## sd(expected1:phase1)                                               0.03 1.00
## sd(expected1:phase2)                                               0.02 1.00
## sd(expected1:difficulty1)                                          0.01 1.00
## sd(expected1:difficulty2)                                          0.02 1.00
## sd(phase1:difficulty1)                                             0.02 1.00
## sd(phase2:difficulty1)                                             0.02 1.00
## sd(phase1:difficulty2)                                             0.02 1.00
## sd(phase2:difficulty2)                                             0.02 1.00
## sd(expected1:phase1:difficulty1)                                   0.02 1.00
## sd(expected1:phase2:difficulty1)                                   0.01 1.00
## sd(expected1:phase1:difficulty2)                                   0.02 1.00
## sd(expected1:phase2:difficulty2)                                   0.02 1.00
## cor(Intercept,expected1)                                           0.50 1.00
## cor(Intercept,phase1)                                              0.31 1.00
## cor(expected1,phase1)                                              0.51 1.00
## cor(Intercept,phase2)                                              0.36 1.00
## cor(expected1,phase2)                                              0.38 1.00
## cor(phase1,phase2)                                                 0.05 1.00
## cor(Intercept,difficulty1)                                         0.10 1.00
## cor(expected1,difficulty1)                                         0.30 1.00
## cor(phase1,difficulty1)                                            0.18 1.00
## cor(phase2,difficulty1)                                            0.30 1.00
## cor(Intercept,difficulty2)                                         0.34 1.00
## cor(expected1,difficulty2)                                         0.30 1.00
## cor(phase1,difficulty2)                                            0.23 1.00
## cor(phase2,difficulty2)                                            0.47 1.00
## cor(difficulty1,difficulty2)                                       0.43 1.00
## cor(Intercept,expected1:phase1)                                    0.36 1.00
## cor(expected1,expected1:phase1)                                    0.53 1.00
## cor(phase1,expected1:phase1)                                       0.44 1.00
## cor(phase2,expected1:phase1)                                       0.21 1.00
## cor(difficulty1,expected1:phase1)                                  0.44 1.00
## cor(difficulty2,expected1:phase1)                                  0.38 1.00
## cor(Intercept,expected1:phase2)                                    0.51 1.00
## cor(expected1,expected1:phase2)                                    0.45 1.00
## cor(phase1,expected1:phase2)                                       0.30 1.00
## cor(phase2,expected1:phase2)                                       0.43 1.00
## cor(difficulty1,expected1:phase2)                                  0.43 1.00
## cor(difficulty2,expected1:phase2)                                  0.52 1.00
## cor(expected1:phase1,expected1:phase2)                             0.41 1.00
## cor(Intercept,expected1:difficulty1)                               0.35 1.00
## cor(expected1,expected1:difficulty1)                               0.45 1.00
## cor(phase1,expected1:difficulty1)                                  0.43 1.00
## cor(phase2,expected1:difficulty1)                                  0.42 1.00
## cor(difficulty1,expected1:difficulty1)                             0.41 1.00
## cor(difficulty2,expected1:difficulty1)                             0.41 1.00
## cor(expected1:phase1,expected1:difficulty1)                        0.45 1.00
## cor(expected1:phase2,expected1:difficulty1)                        0.43 1.00
## cor(Intercept,expected1:difficulty2)                               0.38 1.00
## cor(expected1,expected1:difficulty2)                               0.45 1.00
## cor(phase1,expected1:difficulty2)                                  0.47 1.00
## cor(phase2,expected1:difficulty2)                                  0.46 1.00
## cor(difficulty1,expected1:difficulty2)                             0.42 1.00
## cor(difficulty2,expected1:difficulty2)                             0.41 1.00
## cor(expected1:phase1,expected1:difficulty2)                        0.46 1.00
## cor(expected1:phase2,expected1:difficulty2)                        0.48 1.00
## cor(expected1:difficulty1,expected1:difficulty2)                   0.42 1.00
## cor(Intercept,phase1:difficulty1)                                  0.31 1.00
## cor(expected1,phase1:difficulty1)                                  0.45 1.00
## cor(phase1,phase1:difficulty1)                                     0.53 1.00
## cor(phase2,phase1:difficulty1)                                     0.31 1.00
## cor(difficulty1,phase1:difficulty1)                                0.47 1.00
## cor(difficulty2,phase1:difficulty1)                                0.36 1.00
## cor(expected1:phase1,phase1:difficulty1)                           0.47 1.00
## cor(expected1:phase2,phase1:difficulty1)                           0.38 1.00
## cor(expected1:difficulty1,phase1:difficulty1)                      0.46 1.00
## cor(expected1:difficulty2,phase1:difficulty1)                      0.42 1.00
## cor(Intercept,phase2:difficulty1)                                  0.35 1.00
## cor(expected1,phase2:difficulty1)                                  0.45 1.00
## cor(phase1,phase2:difficulty1)                                     0.31 1.00
## cor(phase2,phase2:difficulty1)                                     0.42 1.00
## cor(difficulty1,phase2:difficulty1)                                0.49 1.00
## cor(difficulty2,phase2:difficulty1)                                0.42 1.00
## cor(expected1:phase1,phase2:difficulty1)                           0.43 1.00
## cor(expected1:phase2,phase2:difficulty1)                           0.48 1.00
## cor(expected1:difficulty1,phase2:difficulty1)                      0.44 1.00
## cor(expected1:difficulty2,phase2:difficulty1)                      0.40 1.00
## cor(phase1:difficulty1,phase2:difficulty1)                         0.40 1.00
## cor(Intercept,phase1:difficulty2)                                  0.45 1.00
## cor(expected1,phase1:difficulty2)                                  0.39 1.00
## cor(phase1,phase1:difficulty2)                                     0.44 1.00
## cor(phase2,phase1:difficulty2)                                     0.36 1.00
## cor(difficulty1,phase1:difficulty2)                                0.39 1.00
## cor(difficulty2,phase1:difficulty2)                                0.43 1.00
## cor(expected1:phase1,phase1:difficulty2)                           0.39 1.00
## cor(expected1:phase2,phase1:difficulty2)                           0.38 1.00
## cor(expected1:difficulty1,phase1:difficulty2)                      0.40 1.00
## cor(expected1:difficulty2,phase1:difficulty2)                      0.42 1.00
## cor(phase1:difficulty1,phase1:difficulty2)                         0.38 1.00
## cor(phase2:difficulty1,phase1:difficulty2)                         0.38 1.00
## cor(Intercept,phase2:difficulty2)                                  0.54 1.00
## cor(expected1,phase2:difficulty2)                                  0.50 1.00
## cor(phase1,phase2:difficulty2)                                     0.28 1.00
## cor(phase2,phase2:difficulty2)                                     0.42 1.00
## cor(difficulty1,phase2:difficulty2)                                0.42 1.00
## cor(difficulty2,phase2:difficulty2)                                0.42 1.00
## cor(expected1:phase1,phase2:difficulty2)                           0.47 1.00
## cor(expected1:phase2,phase2:difficulty2)                           0.56 1.00
## cor(expected1:difficulty1,phase2:difficulty2)                      0.41 1.00
## cor(expected1:difficulty2,phase2:difficulty2)                      0.41 1.00
## cor(phase1:difficulty1,phase2:difficulty2)                         0.35 1.00
## cor(phase2:difficulty1,phase2:difficulty2)                         0.39 1.00
## cor(phase1:difficulty2,phase2:difficulty2)                         0.38 1.00
## cor(Intercept,expected1:phase1:difficulty1)                        0.46 1.00
## cor(expected1,expected1:phase1:difficulty1)                        0.47 1.00
## cor(phase1,expected1:phase1:difficulty1)                           0.54 1.00
## cor(phase2,expected1:phase1:difficulty1)                           0.41 1.00
## cor(difficulty1,expected1:phase1:difficulty1)                      0.41 1.00
## cor(difficulty2,expected1:phase1:difficulty1)                      0.40 1.00
## cor(expected1:phase1,expected1:phase1:difficulty1)                 0.45 1.00
## cor(expected1:phase2,expected1:phase1:difficulty1)                 0.43 1.00
## cor(expected1:difficulty1,expected1:phase1:difficulty1)            0.44 1.00
## cor(expected1:difficulty2,expected1:phase1:difficulty1)            0.46 1.00
## cor(phase1:difficulty1,expected1:phase1:difficulty1)               0.39 1.00
## cor(phase2:difficulty1,expected1:phase1:difficulty1)               0.39 1.00
## cor(phase1:difficulty2,expected1:phase1:difficulty1)               0.41 1.00
## cor(phase2:difficulty2,expected1:phase1:difficulty1)               0.42 1.00
## cor(Intercept,expected1:phase2:difficulty1)                        0.44 1.00
## cor(expected1,expected1:phase2:difficulty1)                        0.50 1.00
## cor(phase1,expected1:phase2:difficulty1)                           0.46 1.00
## cor(phase2,expected1:phase2:difficulty1)                           0.46 1.00
## cor(difficulty1,expected1:phase2:difficulty1)                      0.38 1.00
## cor(difficulty2,expected1:phase2:difficulty1)                      0.39 1.00
## cor(expected1:phase1,expected1:phase2:difficulty1)                 0.41 1.00
## cor(expected1:phase2,expected1:phase2:difficulty1)                 0.41 1.00
## cor(expected1:difficulty1,expected1:phase2:difficulty1)            0.41 1.00
## cor(expected1:difficulty2,expected1:phase2:difficulty1)            0.42 1.00
## cor(phase1:difficulty1,expected1:phase2:difficulty1)               0.40 1.00
## cor(phase2:difficulty1,expected1:phase2:difficulty1)               0.38 1.00
## cor(phase1:difficulty2,expected1:phase2:difficulty1)               0.42 1.00
## cor(phase2:difficulty2,expected1:phase2:difficulty1)               0.41 1.00
## cor(expected1:phase1:difficulty1,expected1:phase2:difficulty1)     0.42 1.00
## cor(Intercept,expected1:phase1:difficulty2)                        0.42 1.00
## cor(expected1,expected1:phase1:difficulty2)                        0.39 1.00
## cor(phase1,expected1:phase1:difficulty2)                           0.42 1.00
## cor(phase2,expected1:phase1:difficulty2)                           0.40 1.00
## cor(difficulty1,expected1:phase1:difficulty2)                      0.39 1.00
## cor(difficulty2,expected1:phase1:difficulty2)                      0.44 1.00
## cor(expected1:phase1,expected1:phase1:difficulty2)                 0.39 1.00
## cor(expected1:phase2,expected1:phase1:difficulty2)                 0.40 1.00
## cor(expected1:difficulty1,expected1:phase1:difficulty2)            0.41 1.00
## cor(expected1:difficulty2,expected1:phase1:difficulty2)            0.43 1.00
## cor(phase1:difficulty1,expected1:phase1:difficulty2)               0.39 1.00
## cor(phase2:difficulty1,expected1:phase1:difficulty2)               0.40 1.00
## cor(phase1:difficulty2,expected1:phase1:difficulty2)               0.41 1.00
## cor(phase2:difficulty2,expected1:phase1:difficulty2)               0.39 1.00
## cor(expected1:phase1:difficulty1,expected1:phase1:difficulty2)     0.41 1.00
## cor(expected1:phase2:difficulty1,expected1:phase1:difficulty2)     0.40 1.00
## cor(Intercept,expected1:phase2:difficulty2)                        0.43 1.00
## cor(expected1,expected1:phase2:difficulty2)                        0.48 1.00
## cor(phase1,expected1:phase2:difficulty2)                           0.48 1.00
## cor(phase2,expected1:phase2:difficulty2)                           0.46 1.00
## cor(difficulty1,expected1:phase2:difficulty2)                      0.39 1.00
## cor(difficulty2,expected1:phase2:difficulty2)                      0.38 1.00
## cor(expected1:phase1,expected1:phase2:difficulty2)                 0.43 1.00
## cor(expected1:phase2,expected1:phase2:difficulty2)                 0.42 1.00
## cor(expected1:difficulty1,expected1:phase2:difficulty2)            0.42 1.00
## cor(expected1:difficulty2,expected1:phase2:difficulty2)            0.42 1.00
## cor(phase1:difficulty1,expected1:phase2:difficulty2)               0.42 1.00
## cor(phase2:difficulty1,expected1:phase2:difficulty2)               0.39 1.00
## cor(phase1:difficulty2,expected1:phase2:difficulty2)               0.40 1.00
## cor(phase2:difficulty2,expected1:phase2:difficulty2)               0.37 1.00
## cor(expected1:phase1:difficulty1,expected1:phase2:difficulty2)     0.44 1.00
## cor(expected1:phase2:difficulty1,expected1:phase2:difficulty2)     0.42 1.00
## cor(expected1:phase1:difficulty2,expected1:phase2:difficulty2)     0.39 1.00
##                                                                Bulk_ESS
## sd(Intercept)                                                      2366
## sd(expected1)                                                      7066
## sd(phase1)                                                         7239
## sd(phase2)                                                         5858
## sd(difficulty1)                                                    4134
## sd(difficulty2)                                                    3095
## sd(expected1:phase1)                                               3029
## sd(expected1:phase2)                                               3068
## sd(expected1:difficulty1)                                          6251
## sd(expected1:difficulty2)                                          3722
## sd(phase1:difficulty1)                                             5189
## sd(phase2:difficulty1)                                             4898
## sd(phase1:difficulty2)                                             6330
## sd(phase2:difficulty2)                                             5160
## sd(expected1:phase1:difficulty1)                                   7427
## sd(expected1:phase2:difficulty1)                                   6694
## sd(expected1:phase1:difficulty2)                                   6824
## sd(expected1:phase2:difficulty2)                                   5243
## cor(Intercept,expected1)                                          14133
## cor(Intercept,phase1)                                              6683
## cor(expected1,phase1)                                              2663
## cor(Intercept,phase2)                                             12625
## cor(expected1,phase2)                                              6010
## cor(phase1,phase2)                                                10967
## cor(Intercept,difficulty1)                                        16432
## cor(expected1,difficulty1)                                        12974
## cor(phase1,difficulty1)                                           15222
## cor(phase2,difficulty1)                                           13511
## cor(Intercept,difficulty2)                                        18051
## cor(expected1,difficulty2)                                        13174
## cor(phase1,difficulty2)                                           14317
## cor(phase2,difficulty2)                                           15036
## cor(difficulty1,difficulty2)                                      15041
## cor(Intercept,expected1:phase1)                                   18603
## cor(expected1,expected1:phase1)                                    9527
## cor(phase1,expected1:phase1)                                      13295
## cor(phase2,expected1:phase1)                                       9736
## cor(difficulty1,expected1:phase1)                                 10090
## cor(difficulty2,expected1:phase1)                                  8434
## cor(Intercept,expected1:phase2)                                   15949
## cor(expected1,expected1:phase2)                                   11900
## cor(phase1,expected1:phase2)                                      14657
## cor(phase2,expected1:phase2)                                      11006
## cor(difficulty1,expected1:phase2)                                  9619
## cor(difficulty2,expected1:phase2)                                  6869
## cor(expected1:phase1,expected1:phase2)                             9188
## cor(Intercept,expected1:difficulty1)                              22552
## cor(expected1,expected1:difficulty1)                              23500
## cor(phase1,expected1:difficulty1)                                 23198
## cor(phase2,expected1:difficulty1)                                 23837
## cor(difficulty1,expected1:difficulty1)                            19936
## cor(difficulty2,expected1:difficulty1)                            17646
## cor(expected1:phase1,expected1:difficulty1)                       16596
## cor(expected1:phase2,expected1:difficulty1)                       18736
## cor(Intercept,expected1:difficulty2)                              20349
## cor(expected1,expected1:difficulty2)                              15884
## cor(phase1,expected1:difficulty2)                                 18050
## cor(phase2,expected1:difficulty2)                                 16314
## cor(difficulty1,expected1:difficulty2)                            15940
## cor(difficulty2,expected1:difficulty2)                            15629
## cor(expected1:phase1,expected1:difficulty2)                       13665
## cor(expected1:phase2,expected1:difficulty2)                       13394
## cor(expected1:difficulty1,expected1:difficulty2)                  12206
## cor(Intercept,phase1:difficulty1)                                 20678
## cor(expected1,phase1:difficulty1)                                 16601
## cor(phase1,phase1:difficulty1)                                    15476
## cor(phase2,phase1:difficulty1)                                    15607
## cor(difficulty1,phase1:difficulty1)                               16586
## cor(difficulty2,phase1:difficulty1)                               15840
## cor(expected1:phase1,phase1:difficulty1)                          17038
## cor(expected1:phase2,phase1:difficulty1)                          17604
## cor(expected1:difficulty1,phase1:difficulty1)                     13039
## cor(expected1:difficulty2,phase1:difficulty1)                     14447
## cor(Intercept,phase2:difficulty1)                                 21087
## cor(expected1,phase2:difficulty1)                                 18351
## cor(phase1,phase2:difficulty1)                                    19609
## cor(phase2,phase2:difficulty1)                                    19276
## cor(difficulty1,phase2:difficulty1)                               14912
## cor(difficulty2,phase2:difficulty1)                               16387
## cor(expected1:phase1,phase2:difficulty1)                          16007
## cor(expected1:phase2,phase2:difficulty1)                          14459
## cor(expected1:difficulty1,phase2:difficulty1)                     13914
## cor(expected1:difficulty2,phase2:difficulty1)                     14616
## cor(phase1:difficulty1,phase2:difficulty1)                        14149
## cor(Intercept,phase1:difficulty2)                                 23002
## cor(expected1,phase1:difficulty2)                                 21002
## cor(phase1,phase1:difficulty2)                                    21507
## cor(phase2,phase1:difficulty2)                                    19294
## cor(difficulty1,phase1:difficulty2)                               18449
## cor(difficulty2,phase1:difficulty2)                               16976
## cor(expected1:phase1,phase1:difficulty2)                          18590
## cor(expected1:phase2,phase1:difficulty2)                          18237
## cor(expected1:difficulty1,phase1:difficulty2)                     12947
## cor(expected1:difficulty2,phase1:difficulty2)                     15489
## cor(phase1:difficulty1,phase1:difficulty2)                        12117
## cor(phase2:difficulty1,phase1:difficulty2)                        12498
## cor(Intercept,phase2:difficulty2)                                 18737
## cor(expected1,phase2:difficulty2)                                 11983
## cor(phase1,phase2:difficulty2)                                    19341
## cor(phase2,phase2:difficulty2)                                    14594
## cor(difficulty1,phase2:difficulty2)                               13320
## cor(difficulty2,phase2:difficulty2)                               11807
## cor(expected1:phase1,phase2:difficulty2)                          12466
## cor(expected1:phase2,phase2:difficulty2)                          10284
## cor(expected1:difficulty1,phase2:difficulty2)                     11650
## cor(expected1:difficulty2,phase2:difficulty2)                     12079
## cor(phase1:difficulty1,phase2:difficulty2)                        11518
## cor(phase2:difficulty1,phase2:difficulty2)                        11069
## cor(phase1:difficulty2,phase2:difficulty2)                        11275
## cor(Intercept,expected1:phase1:difficulty1)                       23070
## cor(expected1,expected1:phase1:difficulty1)                       19927
## cor(phase1,expected1:phase1:difficulty1)                          19826
## cor(phase2,expected1:phase1:difficulty1)                          23301
## cor(difficulty1,expected1:phase1:difficulty1)                     19693
## cor(difficulty2,expected1:phase1:difficulty1)                     18241
## cor(expected1:phase1,expected1:phase1:difficulty1)                20540
## cor(expected1:phase2,expected1:phase1:difficulty1)                18918
## cor(expected1:difficulty1,expected1:phase1:difficulty1)           14901
## cor(expected1:difficulty2,expected1:phase1:difficulty1)           15363
## cor(phase1:difficulty1,expected1:phase1:difficulty1)              13534
## cor(phase2:difficulty1,expected1:phase1:difficulty1)              12602
## cor(phase1:difficulty2,expected1:phase1:difficulty1)              12516
## cor(phase2:difficulty2,expected1:phase1:difficulty1)              15330
## cor(Intercept,expected1:phase2:difficulty1)                       23415
## cor(expected1,expected1:phase2:difficulty1)                       19618
## cor(phase1,expected1:phase2:difficulty1)                          22784
## cor(phase2,expected1:phase2:difficulty1)                          20412
## cor(difficulty1,expected1:phase2:difficulty1)                     19912
## cor(difficulty2,expected1:phase2:difficulty1)                     17690
## cor(expected1:phase1,expected1:phase2:difficulty1)                17898
## cor(expected1:phase2,expected1:phase2:difficulty1)                19036
## cor(expected1:difficulty1,expected1:phase2:difficulty1)           14669
## cor(expected1:difficulty2,expected1:phase2:difficulty1)           16225
## cor(phase1:difficulty1,expected1:phase2:difficulty1)              14340
## cor(phase2:difficulty1,expected1:phase2:difficulty1)              12139
## cor(phase1:difficulty2,expected1:phase2:difficulty1)              12162
## cor(phase2:difficulty2,expected1:phase2:difficulty1)              13742
## cor(expected1:phase1:difficulty1,expected1:phase2:difficulty1)    11303
## cor(Intercept,expected1:phase1:difficulty2)                       24977
## cor(expected1,expected1:phase1:difficulty2)                       21984
## cor(phase1,expected1:phase1:difficulty2)                          25112
## cor(phase2,expected1:phase1:difficulty2)                          23846
## cor(difficulty1,expected1:phase1:difficulty2)                     18818
## cor(difficulty2,expected1:phase1:difficulty2)                     16288
## cor(expected1:phase1,expected1:phase1:difficulty2)                18312
## cor(expected1:phase2,expected1:phase1:difficulty2)                18233
## cor(expected1:difficulty1,expected1:phase1:difficulty2)           13677
## cor(expected1:difficulty2,expected1:phase1:difficulty2)           14578
## cor(phase1:difficulty1,expected1:phase1:difficulty2)              13673
## cor(phase2:difficulty1,expected1:phase1:difficulty2)              13249
## cor(phase1:difficulty2,expected1:phase1:difficulty2)              12764
## cor(phase2:difficulty2,expected1:phase1:difficulty2)              15601
## cor(expected1:phase1:difficulty1,expected1:phase1:difficulty2)    11587
## cor(expected1:phase2:difficulty1,expected1:phase1:difficulty2)    10509
## cor(Intercept,expected1:phase2:difficulty2)                       22721
## cor(expected1,expected1:phase2:difficulty2)                       18641
## cor(phase1,expected1:phase2:difficulty2)                          21970
## cor(phase2,expected1:phase2:difficulty2)                          17872
## cor(difficulty1,expected1:phase2:difficulty2)                     19510
## cor(difficulty2,expected1:phase2:difficulty2)                     16470
## cor(expected1:phase1,expected1:phase2:difficulty2)                19280
## cor(expected1:phase2,expected1:phase2:difficulty2)                17766
## cor(expected1:difficulty1,expected1:phase2:difficulty2)           14471
## cor(expected1:difficulty2,expected1:phase2:difficulty2)           14774
## cor(phase1:difficulty1,expected1:phase2:difficulty2)              15075
## cor(phase2:difficulty1,expected1:phase2:difficulty2)              13598
## cor(phase1:difficulty2,expected1:phase2:difficulty2)              12355
## cor(phase2:difficulty2,expected1:phase2:difficulty2)              13156
## cor(expected1:phase1:difficulty1,expected1:phase2:difficulty2)    11294
## cor(expected1:phase2:difficulty1,expected1:phase2:difficulty2)    11446
## cor(expected1:phase1:difficulty2,expected1:phase2:difficulty2)    10626
##                                                                Tail_ESS
## sd(Intercept)                                                      4846
## sd(expected1)                                                      7666
## sd(phase1)                                                        11628
## sd(phase2)                                                         7273
## sd(difficulty1)                                                    2776
## sd(difficulty2)                                                    3998
## sd(expected1:phase1)                                               2182
## sd(expected1:phase2)                                               2836
## sd(expected1:difficulty1)                                          7980
## sd(expected1:difficulty2)                                          5876
## sd(phase1:difficulty1)                                             6095
## sd(phase2:difficulty1)                                             6781
## sd(phase1:difficulty2)                                             6790
## sd(phase2:difficulty2)                                             5173
## sd(expected1:phase1:difficulty1)                                   6748
## sd(expected1:phase2:difficulty1)                                   7628
## sd(expected1:phase1:difficulty2)                                   7751
## sd(expected1:phase2:difficulty2)                                   7346
## cor(Intercept,expected1)                                          13901
## cor(Intercept,phase1)                                             10052
## cor(expected1,phase1)                                              5453
## cor(Intercept,phase2)                                             13120
## cor(expected1,phase2)                                              8180
## cor(phase1,phase2)                                                13010
## cor(Intercept,difficulty1)                                        10045
## cor(expected1,difficulty1)                                        13322
## cor(phase1,difficulty1)                                           11806
## cor(phase2,difficulty1)                                           12642
## cor(Intercept,difficulty2)                                        12744
## cor(expected1,difficulty2)                                        13916
## cor(phase1,difficulty2)                                           11596
## cor(phase2,difficulty2)                                           13506
## cor(difficulty1,difficulty2)                                      13240
## cor(Intercept,expected1:phase1)                                   13281
## cor(expected1,expected1:phase1)                                   11700
## cor(phase1,expected1:phase1)                                      12199
## cor(phase2,expected1:phase1)                                      12181
## cor(difficulty1,expected1:phase1)                                 11920
## cor(difficulty2,expected1:phase1)                                 13198
## cor(Intercept,expected1:phase2)                                   12123
## cor(expected1,expected1:phase2)                                   13146
## cor(phase1,expected1:phase2)                                      13417
## cor(phase2,expected1:phase2)                                      12579
## cor(difficulty1,expected1:phase2)                                 11882
## cor(difficulty2,expected1:phase2)                                 11126
## cor(expected1:phase1,expected1:phase2)                            13013
## cor(Intercept,expected1:difficulty1)                              13046
## cor(expected1,expected1:difficulty1)                              13311
## cor(phase1,expected1:difficulty1)                                 14136
## cor(phase2,expected1:difficulty1)                                 13834
## cor(difficulty1,expected1:difficulty1)                            13714
## cor(difficulty2,expected1:difficulty1)                            14953
## cor(expected1:phase1,expected1:difficulty1)                       14001
## cor(expected1:phase2,expected1:difficulty1)                       14017
## cor(Intercept,expected1:difficulty2)                              13541
## cor(expected1,expected1:difficulty2)                              13431
## cor(phase1,expected1:difficulty2)                                 13038
## cor(phase2,expected1:difficulty2)                                 13622
## cor(difficulty1,expected1:difficulty2)                            12615
## cor(difficulty2,expected1:difficulty2)                            13899
## cor(expected1:phase1,expected1:difficulty2)                       14302
## cor(expected1:phase2,expected1:difficulty2)                       14382
## cor(expected1:difficulty1,expected1:difficulty2)                  14845
## cor(Intercept,phase1:difficulty1)                                 13282
## cor(expected1,phase1:difficulty1)                                 13648
## cor(phase1,phase1:difficulty1)                                    12386
## cor(phase2,phase1:difficulty1)                                    13542
## cor(difficulty1,phase1:difficulty1)                               13964
## cor(difficulty2,phase1:difficulty1)                               14759
## cor(expected1:phase1,phase1:difficulty1)                          14109
## cor(expected1:phase2,phase1:difficulty1)                          13794
## cor(expected1:difficulty1,phase1:difficulty1)                     14056
## cor(expected1:difficulty2,phase1:difficulty1)                     13721
## cor(Intercept,phase2:difficulty1)                                 13027
## cor(expected1,phase2:difficulty1)                                 13619
## cor(phase1,phase2:difficulty1)                                    13105
## cor(phase2,phase2:difficulty1)                                    13993
## cor(difficulty1,phase2:difficulty1)                               13596
## cor(difficulty2,phase2:difficulty1)                               14732
## cor(expected1:phase1,phase2:difficulty1)                          14957
## cor(expected1:phase2,phase2:difficulty1)                          14145
## cor(expected1:difficulty1,phase2:difficulty1)                     13741
## cor(expected1:difficulty2,phase2:difficulty1)                     15255
## cor(phase1:difficulty1,phase2:difficulty1)                        15542
## cor(Intercept,phase1:difficulty2)                                 13489
## cor(expected1,phase1:difficulty2)                                 13678
## cor(phase1,phase1:difficulty2)                                    12851
## cor(phase2,phase1:difficulty2)                                    13674
## cor(difficulty1,phase1:difficulty2)                               13341
## cor(difficulty2,phase1:difficulty2)                               13816
## cor(expected1:phase1,phase1:difficulty2)                          13784
## cor(expected1:phase2,phase1:difficulty2)                          14706
## cor(expected1:difficulty1,phase1:difficulty2)                     13261
## cor(expected1:difficulty2,phase1:difficulty2)                     14463
## cor(phase1:difficulty1,phase1:difficulty2)                        14085
## cor(phase2:difficulty1,phase1:difficulty2)                        14235
## cor(Intercept,phase2:difficulty2)                                 12832
## cor(expected1,phase2:difficulty2)                                 13495
## cor(phase1,phase2:difficulty2)                                    13723
## cor(phase2,phase2:difficulty2)                                    13388
## cor(difficulty1,phase2:difficulty2)                               13522
## cor(difficulty2,phase2:difficulty2)                               13222
## cor(expected1:phase1,phase2:difficulty2)                          13607
## cor(expected1:phase2,phase2:difficulty2)                          12021
## cor(expected1:difficulty1,phase2:difficulty2)                     14430
## cor(expected1:difficulty2,phase2:difficulty2)                     13639
## cor(phase1:difficulty1,phase2:difficulty2)                        13764
## cor(phase2:difficulty1,phase2:difficulty2)                        13464
## cor(phase1:difficulty2,phase2:difficulty2)                        15136
## cor(Intercept,expected1:phase1:difficulty1)                       13700
## cor(expected1,expected1:phase1:difficulty1)                       13061
## cor(phase1,expected1:phase1:difficulty1)                          13947
## cor(phase2,expected1:phase1:difficulty1)                          13971
## cor(difficulty1,expected1:phase1:difficulty1)                     14156
## cor(difficulty2,expected1:phase1:difficulty1)                     13666
## cor(expected1:phase1,expected1:phase1:difficulty1)                13708
## cor(expected1:phase2,expected1:phase1:difficulty1)                14229
## cor(expected1:difficulty1,expected1:phase1:difficulty1)           14680
## cor(expected1:difficulty2,expected1:phase1:difficulty1)           14001
## cor(phase1:difficulty1,expected1:phase1:difficulty1)              13316
## cor(phase2:difficulty1,expected1:phase1:difficulty1)              14134
## cor(phase1:difficulty2,expected1:phase1:difficulty1)              14264
## cor(phase2:difficulty2,expected1:phase1:difficulty1)              15749
## cor(Intercept,expected1:phase2:difficulty1)                       13120
## cor(expected1,expected1:phase2:difficulty1)                       13801
## cor(phase1,expected1:phase2:difficulty1)                          13974
## cor(phase2,expected1:phase2:difficulty1)                          14136
## cor(difficulty1,expected1:phase2:difficulty1)                     14253
## cor(difficulty2,expected1:phase2:difficulty1)                     14454
## cor(expected1:phase1,expected1:phase2:difficulty1)                13891
## cor(expected1:phase2,expected1:phase2:difficulty1)                14646
## cor(expected1:difficulty1,expected1:phase2:difficulty1)           14325
## cor(expected1:difficulty2,expected1:phase2:difficulty1)           15700
## cor(phase1:difficulty1,expected1:phase2:difficulty1)              14899
## cor(phase2:difficulty1,expected1:phase2:difficulty1)              14768
## cor(phase1:difficulty2,expected1:phase2:difficulty1)              13772
## cor(phase2:difficulty2,expected1:phase2:difficulty1)              14138
## cor(expected1:phase1:difficulty1,expected1:phase2:difficulty1)    14170
## cor(Intercept,expected1:phase1:difficulty2)                       13094
## cor(expected1,expected1:phase1:difficulty2)                       13509
## cor(phase1,expected1:phase1:difficulty2)                          12429
## cor(phase2,expected1:phase1:difficulty2)                          14026
## cor(difficulty1,expected1:phase1:difficulty2)                     14670
## cor(difficulty2,expected1:phase1:difficulty2)                     14558
## cor(expected1:phase1,expected1:phase1:difficulty2)                14823
## cor(expected1:phase2,expected1:phase1:difficulty2)                14127
## cor(expected1:difficulty1,expected1:phase1:difficulty2)           14486
## cor(expected1:difficulty2,expected1:phase1:difficulty2)           14711
## cor(phase1:difficulty1,expected1:phase1:difficulty2)              14330
## cor(phase2:difficulty1,expected1:phase1:difficulty2)              14702
## cor(phase1:difficulty2,expected1:phase1:difficulty2)              15253
## cor(phase2:difficulty2,expected1:phase1:difficulty2)              15659
## cor(expected1:phase1:difficulty1,expected1:phase1:difficulty2)    14662
## cor(expected1:phase2:difficulty1,expected1:phase1:difficulty2)    13875
## cor(Intercept,expected1:phase2:difficulty2)                       13616
## cor(expected1,expected1:phase2:difficulty2)                       14173
## cor(phase1,expected1:phase2:difficulty2)                          13712
## cor(phase2,expected1:phase2:difficulty2)                          13197
## cor(difficulty1,expected1:phase2:difficulty2)                     14104
## cor(difficulty2,expected1:phase2:difficulty2)                     14255
## cor(expected1:phase1,expected1:phase2:difficulty2)                14359
## cor(expected1:phase2,expected1:phase2:difficulty2)                14830
## cor(expected1:difficulty1,expected1:phase2:difficulty2)           15253
## cor(expected1:difficulty2,expected1:phase2:difficulty2)           14510
## cor(phase1:difficulty1,expected1:phase2:difficulty2)              14737
## cor(phase2:difficulty1,expected1:phase2:difficulty2)              14498
## cor(phase1:difficulty2,expected1:phase2:difficulty2)              13653
## cor(phase2:difficulty2,expected1:phase2:difficulty2)              14074
## cor(expected1:phase1:difficulty1,expected1:phase2:difficulty2)    14532
## cor(expected1:phase2:difficulty1,expected1:phase2:difficulty2)    14123
## cor(expected1:phase1:difficulty2,expected1:phase2:difficulty2)    14521
## 
## ~trl (Number of levels: 288) 
##                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(Intercept)                  0.07      0.00     0.06     0.08 1.00     6048
## sd(diagnosis1)                 0.01      0.00     0.00     0.02 1.00     3418
## sd(diagnosis2)                 0.00      0.00     0.00     0.01 1.00     4559
## cor(Intercept,diagnosis1)      0.49      0.23    -0.03     0.87 1.00    10682
## cor(Intercept,diagnosis2)     -0.03      0.37    -0.72     0.69 1.00    18822
## cor(diagnosis1,diagnosis2)    -0.23      0.42    -0.88     0.65 1.00     8090
##                            Tail_ESS
## sd(Intercept)                 10131
## sd(diagnosis1)                 3325
## sd(diagnosis2)                 7185
## cor(Intercept,diagnosis1)      9190
## cor(Intercept,diagnosis2)     12435
## cor(diagnosis1,diagnosis2)    12092
## 
## Regression Coefficients:
##                                         Estimate Est.Error l-95% CI u-95% CI
## Intercept                                   6.17      0.03     6.12     6.22
## diagnosis1                                  0.02      0.03    -0.03     0.07
## diagnosis2                                  0.02      0.03    -0.03     0.07
## expected1                                  -0.04      0.01    -0.05    -0.03
## phase1                                      0.02      0.01    -0.01     0.04
## phase2                                     -0.01      0.01    -0.02     0.01
## difficulty1                                -0.03      0.01    -0.05    -0.02
## difficulty2                                -0.00      0.01    -0.02     0.01
## diagnosis1:expected1                       -0.01      0.01    -0.02     0.00
## diagnosis2:expected1                       -0.00      0.00    -0.01     0.01
## diagnosis1:phase1                          -0.01      0.01    -0.03     0.01
## diagnosis2:phase1                          -0.00      0.01    -0.02     0.02
## diagnosis1:phase2                          -0.01      0.01    -0.02     0.01
## diagnosis2:phase2                           0.00      0.01    -0.01     0.02
## expected1:phase1                           -0.04      0.01    -0.06    -0.02
## expected1:phase2                            0.00      0.01    -0.01     0.02
## diagnosis1:difficulty1                     -0.01      0.01    -0.02     0.00
## diagnosis2:difficulty1                      0.00      0.01    -0.01     0.02
## diagnosis1:difficulty2                      0.02      0.01     0.01     0.03
## diagnosis2:difficulty2                     -0.02      0.01    -0.03    -0.01
## expected1:difficulty1                      -0.01      0.01    -0.03     0.00
## expected1:difficulty2                       0.00      0.01    -0.01     0.02
## phase1:difficulty1                         -0.01      0.01    -0.03     0.02
## phase2:difficulty1                         -0.00      0.01    -0.02     0.02
## phase1:difficulty2                          0.01      0.01    -0.01     0.04
## phase2:difficulty2                         -0.01      0.01    -0.03     0.01
## diagnosis1:expected1:phase1                -0.00      0.01    -0.02     0.01
## diagnosis2:expected1:phase1                 0.01      0.01    -0.00     0.02
## diagnosis1:expected1:phase2                 0.00      0.01    -0.01     0.01
## diagnosis2:expected1:phase2                -0.00      0.01    -0.01     0.01
## diagnosis1:expected1:difficulty1           -0.01      0.01    -0.02     0.01
## diagnosis2:expected1:difficulty1            0.01      0.01    -0.00     0.02
## diagnosis1:expected1:difficulty2           -0.00      0.01    -0.01     0.01
## diagnosis2:expected1:difficulty2            0.00      0.01    -0.01     0.01
## diagnosis1:phase1:difficulty1               0.00      0.01    -0.01     0.02
## diagnosis2:phase1:difficulty1               0.00      0.01    -0.01     0.02
## diagnosis1:phase2:difficulty1              -0.01      0.01    -0.02     0.00
## diagnosis2:phase2:difficulty1              -0.00      0.01    -0.02     0.01
## diagnosis1:phase1:difficulty2               0.00      0.01    -0.01     0.02
## diagnosis2:phase1:difficulty2              -0.01      0.01    -0.03     0.00
## diagnosis1:phase2:difficulty2              -0.00      0.01    -0.02     0.01
## diagnosis2:phase2:difficulty2               0.01      0.01    -0.00     0.03
## expected1:phase1:difficulty1               -0.01      0.01    -0.03     0.02
## expected1:phase2:difficulty1               -0.00      0.01    -0.02     0.02
## expected1:phase1:difficulty2                0.01      0.01    -0.01     0.03
## expected1:phase2:difficulty2               -0.01      0.01    -0.03     0.01
## diagnosis1:expected1:phase1:difficulty1     0.00      0.01    -0.01     0.02
## diagnosis2:expected1:phase1:difficulty1     0.00      0.01    -0.01     0.02
## diagnosis1:expected1:phase2:difficulty1     0.01      0.01    -0.00     0.02
## diagnosis2:expected1:phase2:difficulty1    -0.01      0.01    -0.02     0.00
## diagnosis1:expected1:phase1:difficulty2     0.00      0.01    -0.01     0.02
## diagnosis2:expected1:phase1:difficulty2    -0.00      0.01    -0.02     0.01
## diagnosis1:expected1:phase2:difficulty2    -0.01      0.01    -0.02     0.00
## diagnosis2:expected1:phase2:difficulty2     0.00      0.01    -0.01     0.01
##                                         Rhat Bulk_ESS Tail_ESS
## Intercept                               1.00     1257     2922
## diagnosis1                              1.00     1384     3130
## diagnosis2                              1.00     1632     3199
## expected1                               1.00     5728     8644
## phase1                                  1.00     5107     8540
## phase2                                  1.00     5216     9337
## difficulty1                             1.00     4814     8336
## difficulty2                             1.00     4858     8799
## diagnosis1:expected1                    1.00    10348    13271
## diagnosis2:expected1                    1.00    10286    12734
## diagnosis1:phase1                       1.00     5440     9226
## diagnosis2:phase1                       1.00     5033     8789
## diagnosis1:phase2                       1.00     9242    11314
## diagnosis2:phase2                       1.00     9484    11579
## expected1:phase1                        1.00     5887     9381
## expected1:phase2                        1.00     5239     9521
## diagnosis1:difficulty1                  1.00     9997    11062
## diagnosis2:difficulty1                  1.00    10562    12419
## diagnosis1:difficulty2                  1.00    11211    12983
## diagnosis2:difficulty2                  1.00    11525    13218
## expected1:difficulty1                   1.00     4894     8458
## expected1:difficulty2                   1.00     5300     8739
## phase1:difficulty1                      1.00     4742     8273
## phase2:difficulty1                      1.00     4199     8491
## phase1:difficulty2                      1.00     4816     9280
## phase2:difficulty2                      1.00     4547     7937
## diagnosis1:expected1:phase1             1.00    13718    14036
## diagnosis2:expected1:phase1             1.00    13722    13021
## diagnosis1:expected1:phase2             1.00    10816    12500
## diagnosis2:expected1:phase2             1.00    11218    13346
## diagnosis1:expected1:difficulty1        1.00    10306    13500
## diagnosis2:expected1:difficulty1        1.00    11979    13393
## diagnosis1:expected1:difficulty2        1.00    11054    12978
## diagnosis2:expected1:difficulty2        1.00    11268    12715
## diagnosis1:phase1:difficulty1           1.00     9380    12504
## diagnosis2:phase1:difficulty1           1.00    10149    12743
## diagnosis1:phase2:difficulty1           1.00     9314    11792
## diagnosis2:phase2:difficulty1           1.00     9906    12839
## diagnosis1:phase1:difficulty2           1.00     9038    12396
## diagnosis2:phase1:difficulty2           1.00     9988    13029
## diagnosis1:phase2:difficulty2           1.00     9117    13078
## diagnosis2:phase2:difficulty2           1.00     9025    12760
## expected1:phase1:difficulty1            1.00     5180     8462
## expected1:phase2:difficulty1            1.00     4292     7684
## expected1:phase1:difficulty2            1.00     5284     8058
## expected1:phase2:difficulty2            1.00     4688     7872
## diagnosis1:expected1:phase1:difficulty1 1.00     9765    13094
## diagnosis2:expected1:phase1:difficulty1 1.00    10073    12791
## diagnosis1:expected1:phase2:difficulty1 1.00     9494    12343
## diagnosis2:expected1:phase2:difficulty1 1.00     9763    12819
## diagnosis1:expected1:phase1:difficulty2 1.00     9227    12066
## diagnosis2:expected1:phase1:difficulty2 1.00     9940    12701
## diagnosis1:expected1:phase2:difficulty2 1.00     9034    13206
## diagnosis2:expected1:phase2:difficulty2 1.00     9861    13223
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.25      0.00     0.24     0.25 1.00    11159    11907
## ndt     106.89      5.16    96.49   116.71 1.00    10792    12313
## 
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# get the estimates and compute group comparisons
df.m.pal = post.draws %>% 
  select(starts_with("b_")) %>%
  mutate(
    b_postvolatile = - b_phase1 - b_phase2,
    b_difficult    = - b_difficulty1 - b_difficulty2
  )

# plot the posterior distributions
df.m.pal %>%
  pivot_longer(cols = starts_with("b_"), names_to = "coef", values_to = "estimate") %>%
  subset(!startsWith(coef, "b_Int")) %>%
  mutate(
    coef = substr(coef, 3, nchar(coef)),
    coef = str_replace_all(coef, ":", " x "),
    coef = str_replace_all(coef, "diagnosis1", "ADHD"),
    coef = str_replace_all(coef, "diagnosis2", "ADHD+ASD"),
    coef = str_replace_all(coef, "expected1", "expected"),
    coef = str_replace_all(coef, "expected2", "unexpected"),
    coef = str_replace_all(coef, "phase1", "prevolatile"),
    coef = str_replace_all(coef, "phase2", "volatile"),
    coef = str_replace_all(coef, "difficulty1", "easy"),
    coef = str_replace_all(coef, "difficulty2", "medium"),
    coef = fct_reorder(coef, desc(estimate))
  ) %>% 
  group_by(coef) %>%
  mutate(
    cred = case_when(
      (mean(estimate) < 0 & quantile(estimate, probs = 0.975) < 0) |
        (mean(estimate) > 0 & quantile(estimate, probs = 0.025) > 0) ~ "credible",
      T ~ "not credible"
    )
  ) %>% ungroup() %>%
  ggplot(aes(x = estimate, y = coef, fill = cred)) +
  geom_vline(xintercept = 0, linetype = 'dashed') +
  ggdist::stat_halfeye(alpha = 0.7) + ylab(NULL) + theme_bw() +
  scale_fill_manual(values = c(c_dark, c_light)) + theme(legend.position = "none")

## create design matrix to figure out how to set comparisons for hypotheses
df.des = cbind(df.pal %>% select(diagnosis, expected, phase, difficulty), 
               model.matrix(~ diagnosis * expected * phase * difficulty, 
                            data = df.pal)) %>% distinct()

# H1c: COMP(unexp-exp) != ADHD(unexp-exp)
as.data.frame(t(df.des %>% 
    filter(diagnosis != "BOTH") %>%
    group_by(diagnosis, expected) %>%
    summarise(across(where(is.numeric), ~ mean(.x))) %>%
    group_by(diagnosis) %>%
    summarise(across(where(is.numeric), ~ diff(.x))) %>% # unexpected - expected
    select(where(is.numeric)) %>%
    map_df(~ diff(.x)))) %>%
  filter(V1 != 0) # COMP - ADHD
##                      V1
## diagnosis1:expected1  4
## diagnosis2:expected1  2
h1c = hypothesis(m.pal, "0 > 2*diagnosis1:expected1 + diagnosis2:expected1")
h1c
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (0)-(2*diagnosis1... > 0     0.01      0.01        0     0.03       9.44
##   Post.Prob Star
## 1       0.9     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
# COMP(unexp-exp) != BOTH(unexp-exp)
e.both = hypothesis(m.pal, "0 > diagnosis1:expected1 + 2*diagnosis2:expected1", alpha = 0.025)
e.both
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (0)-(diagnosis1:e... > 0     0.01      0.01    -0.01     0.02       3.49
##   Post.Prob Star
## 1      0.78     
## ---
## 'CI': 95%-CI for one-sided and 97.5%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 97.5%;
## for two-sided hypotheses, the value tested against lies outside the 97.5%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
# ADHD(unexp-exp) != BOTH(unexp-exp)
e.clin = hypothesis(m.pal, "diagnosis1:expected1 < diagnosis2:expected1", alpha = 0.025)
e.clin
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (diagnosis1:expec... < 0        0      0.01    -0.02     0.01        2.5
##   Post.Prob Star
## 1      0.71     
## ---
## 'CI': 95%-CI for one-sided and 97.5%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 97.5%;
## for two-sided hypotheses, the value tested against lies outside the 97.5%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
## exploration: task effects

# expected versus unexpected
e1 = hypothesis(m.pal, "0 > 2*expected1", alpha = 0.025)
e1
## Hypothesis Tests for class b:
##              Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (0)-(2*expected1) > 0     0.08      0.01     0.06     0.11        Inf
##   Post.Prob Star
## 1         1    *
## ---
## 'CI': 95%-CI for one-sided and 97.5%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 97.5%;
## for two-sided hypotheses, the value tested against lies outside the 97.5%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
# volatile versus prevolatile
e2 = hypothesis(m.pal, "0 < phase1 - phase2", alpha = 0.025)
e2
## Hypothesis Tests for class b:
##                Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (0)-(phase1-phase2) < 0    -0.02      0.02    -0.06     0.01      15.03
##   Post.Prob Star
## 1      0.94     
## ---
## 'CI': 95%-CI for one-sided and 97.5%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 97.5%;
## for two-sided hypotheses, the value tested against lies outside the 97.5%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
# volatile versus postvolatile
e3 = hypothesis(m.pal, "0 < phase1 + 2*phase2", alpha = 0.025)
e3
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (0)-(phase1+2*pha... < 0        0      0.02    -0.03     0.03       1.15
##   Post.Prob Star
## 1      0.54     
## ---
## 'CI': 95%-CI for one-sided and 97.5%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 97.5%;
## for two-sided hypotheses, the value tested against lies outside the 97.5%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
# difficult versus medium
e4 = hypothesis(m.pal, "0 < -difficulty1 - 2*difficulty2", alpha = 0.025)
e4
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (0)-(-difficulty1... < 0    -0.04      0.02    -0.07    -0.01     366.35
##   Post.Prob Star
## 1         1    *
## ---
## 'CI': 95%-CI for one-sided and 97.5%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 97.5%;
## for two-sided hypotheses, the value tested against lies outside the 97.5%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
# medium versus easy
e5 = hypothesis(m.pal, "0 < -difficulty1 + difficulty2", alpha = 0.025)
e5
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (0)-(-difficulty1... < 0    -0.03      0.01    -0.06        0      44.92
##   Post.Prob Star
## 1      0.98    *
## ---
## 'CI': 95%-CI for one-sided and 97.5%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 97.5%;
## for two-sided hypotheses, the value tested against lies outside the 97.5%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
# difficult versus easy
e6 = hypothesis(m.pal, "0 < -2*difficulty1 - difficulty2", alpha = 0.025)
e6
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (0)-(-2*difficult... < 0    -0.07      0.02     -0.1    -0.04        Inf
##   Post.Prob Star
## 1         1    *
## ---
## 'CI': 95%-CI for one-sided and 97.5%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 97.5%;
## for two-sided hypotheses, the value tested against lies outside the 97.5%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
## extract predicted differences in ms instead of log data
df.new = df.pal %>% 
  select(diagnosis, phase, expected, difficulty) %>% 
  distinct() %>%
  mutate(
    condition = paste(diagnosis, phase, expected, difficulty, sep = "_")
  )
df.ms = as.data.frame(
  fitted(m.pal, summary = F, 
               newdata = df.new %>% select(diagnosis, phase, expected, difficulty), 
               re_formula = NA))
colnames(df.ms) = df.new$condition

# calculate our difference columns
df.ms = df.ms %>%
  mutate(
    COMP_expected     = rowMeans(across(matches("COMP_.*_expected_.*"))),
    COMP_unexpected   = rowMeans(across(matches("COMP_.*_unexpected_.*"))),
    ADHD_expected      = rowMeans(across(matches("ADHD_.*_expected_.*"))),
    ADHD_unexpected    = rowMeans(across(matches("ADHD_.*_unexpected_.*"))),
    BOTH_expected      = rowMeans(across(matches("BOTH_.*_expected_.*"))),
    BOTH_unexpected    = rowMeans(across(matches("BOTH_.*_unexpected_.*"))),
    COMP_unexp_exp    = COMP_unexpected - COMP_expected,
    ADHD_unexp_exp     = ADHD_unexpected  - ADHD_expected,
    BOTH_unexp_exp     = BOTH_unexpected  - BOTH_expected,
    `h1c_COMP-ADHD unexpected-expected`  = ADHD_unexp_exp - COMP_unexp_exp,
    `e_COMP-BOTH unexpected-expected`  = BOTH_unexp_exp - COMP_unexp_exp,
    `e_ADHD-BOTH unexpected-expected`  = ADHD_unexp_exp - BOTH_unexp_exp,
    `e1_unexpected-expected`           = rowMeans(across(matches(".*_unexpected_.*"))) - 
      rowMeans(across(matches(".*_expected_.*"))),
    `e2_volatile-prevolatile`       = rowMeans(across(matches(".*_volatile_.*"))) - 
      rowMeans(across(matches(".*_prevolatile_.*"))),
    `e3_volatile-postvolatile`      = rowMeans(across(matches(".*_volatile_.*"))) - 
      rowMeans(across(matches(".*_postvolatile_.*"))),
    `e4_difficult-medium`      = rowMeans(across(matches(".*_difficult"))) - 
      rowMeans(across(matches(".*_medium"))),
    `e5_medium-easy`      = rowMeans(across(matches(".*_medium"))) - 
      rowMeans(across(matches(".*_easy"))),
    `e6_difficult-easy`     = rowMeans(across(matches(".*_difficult"))) - 
      rowMeans(across(matches(".*_easy")))
  )

kable(df.ms %>% 
  select(starts_with("e") | starts_with("h")) %>% 
  summarise_all(list(lower.ci = lower_ci, mean = mean, upper.ci = upper_ci)) %>%
  t %>% 
  as.data.frame %>% 
  rownames_to_column() %>%
  separate(rowname, into = c("id", "comparison", "stat"), sep = "_") %>%
  pivot_wider(names_from = stat, values_from = V1) %>%
  arrange(comparison), digits = 2)
id comparison lower.ci mean upper.ci
e ADHD-BOTH unexpected-expected -12.73 4.50 21.66
h1c COMP-ADHD unexpected-expected -4.37 12.89 30.22
e COMP-BOTH unexpected-expected -8.61 8.39 25.57
e6 difficult-easy -11.55 3.47 18.97
e4 difficult-medium 5.59 20.76 36.26
e5 medium-easy -0.26 14.23 28.62
e1 unexpected-expected 27.75 40.71 54.22
e3 volatile-postvolatile -14.98 0.80 16.60
e2 volatile-prevolatile -29.65 -13.41 2.76
equivalence_test(df.ms %>% select(`h1c_COMP-ADHD unexpected-expected`,
                                  `e_COMP-BOTH unexpected-expected`,
                                  `e_ADHD-BOTH unexpected-expected`),
                 range = rope_range(m.pal))
## # Test for Practical Equivalence
## 
##   ROPE: [-17.71 17.71]
## 
## Parameter                         |        H0 | inside ROPE |         95% HDI
## -----------------------------------------------------------------------------
## h1c_COMP-ADHD unexpected-expected | Undecided |     71.64 % |  [-4.37, 30.22]
## e_COMP-BOTH unexpected-expected   | Undecided |     88.00 % |  [-8.61, 25.57]
## e_ADHD-BOTH unexpected-expected   | Undecided |     95.87 % | [-12.73, 21.66]
# calculate effect sizes
df.effect = post.draws %>%
  mutate(across(starts_with("sd")|starts_with("sigma"), ~.^2)) %>%
  mutate(
    sumvar = sqrt(rowSums(select(., starts_with("sd")|starts_with("sigma")))),
    h1c    = -(2*`b_diagnosis1:expected1` + `b_diagnosis2:expected1`) / sumvar,
    e1     = -(2*b_expected1) / sumvar,
    e6     = -(2*`b_difficulty1` + `b_difficulty2`) / sumvar
  )

kable(df.effect %>% select(starts_with("e")|starts_with("h")) %>%
        pivot_longer(cols = everything(), values_to = "estimate") %>%
        group_by(name) %>%
        summarise(
          ci.lo = lower_ci(estimate),
          mean  = mean(estimate),
          ci.hi = upper_ci(estimate),
          interpret = interpret_cohens_d(mean)
        ), digits = 3
)
name ci.lo mean ci.hi interpret
e1 0.169 0.249 0.331 small
e6 0.130 0.219 0.312 small
h1c -0.018 0.035 0.089 very small

h1c: estimate = 0.01 [0, 0.03], posterior probability = 90.42%

3.3 Plots

# rain cloud plot including all factors

df.pal %>%
  ggplot(aes(diagnosis, rt.cor, fill = expected, colour = expected)) + #
  geom_rain(rain.side = 'r',
boxplot.args = list(color = "black", outlier.shape = NA, show.legend = FALSE, alpha = .8),
violin.args = list(color = "black", outlier.shape = NA, alpha = .8),
boxplot.args.pos = list(
  position = ggpp::position_dodgenudge(x = 0, width = 0.3), width = 0.3
),
point.args = list(show_guide = FALSE, alpha = .5),
violin.args.pos = list(
  width = 0.6, position = position_nudge(x = 0.16)),
point.args.pos = list(position = ggpp::position_dodgenudge(x = -0.25, width = 0.1))) +
  scale_fill_manual(values = col.cnd) +
  scale_color_manual(values = col.cnd) +
  #ylim(0, 1) +
  labs(title = "Reaction times per subject", x = "", y = "rt (ms)") +
  facet_wrap(difficulty ~ phase) +
  theme_bw() + 
  theme(legend.position = "bottom", plot.title = element_text(hjust = 0.5), 
        legend.direction = "horizontal", text = element_text(size = 15))

# rain cloud plot excluding difficulty
df.pal %>%
  group_by(subID, diagnosis, expected, phase) %>%
  summarise(
    rt.cor = median(rt.cor, na.rm = T)
  ) %>%
  group_by(subID, diagnosis, phase) %>%
  summarise(
    rt.exp = diff(rt.cor)
  ) %>%
  ggplot(aes(phase, rt.exp, fill = diagnosis, colour = diagnosis)) + #
  geom_rain(rain.side = 'r',
boxplot.args = list(color = "black", outlier.shape = NA, show.legend = FALSE, alpha = .8),
violin.args = list(color = "black", outlier.shape = NA, alpha = .8),
boxplot.args.pos = list(
  position = ggpp::position_dodgenudge(x = 0, width = 0.3), width = 0.3
),
point.args = list(show_guide = FALSE, alpha = .5),
violin.args.pos = list(
  width = 0.6, position = position_nudge(x = 0.16)),
point.args.pos = list(position = ggpp::position_dodgenudge(x = -0.25, width = 0.1))) +
  scale_fill_manual(values = col.grp) +
  scale_color_manual(values = col.grp) +
  #ylim(0, 1) +
  labs(title = "", x = "", y = "rt difference (ms)") +
  theme_bw() + 
  theme(legend.position = "bottom", plot.title = element_blank(), 
        legend.direction = "horizontal", text = element_text(size = 15),
        legend.title = element_blank())

ggsave("plots/FigRT.svg", units = "cm", width = 27, height = 9)

# plot transition effects specifically
df.pal %>% 
  group_by(subID, diagnosis, phase, expected) %>%
  summarise(rt.cor = median(rt.cor, na.rm = T)) %>%
  group_by(diagnosis, phase, expected) %>%
  summarise(
    rt.mn = mean(rt.cor),
    rt.se = sd(rt.cor)
  ) %>%
  ggplot(aes(y = rt.mn, x = phase, 
             group = diagnosis, colour = diagnosis)) +
  geom_line(position = position_dodge(0.4), linewidth = 1) +
  geom_errorbar(aes(ymin = rt.mn - rt.se, 
                    ymax = rt.mn + rt.se), linewidth = 1, width = 0.5, 
                position = position_dodge(0.4)) + 
  geom_point(position = position_dodge(0.4), size = 5) + 
  labs (x = "phase", y = "rt (ms)") + 
  ggtitle("Diagnosis x phase x expected") + 
  facet_wrap(. ~ expected) +
  scale_fill_manual(values = col.grp) +
  scale_color_manual(values = col.grp) +
  theme_bw() + 
  theme(legend.position = "bottom", plot.title = element_text(hjust = 0.5))

4 Explorative analysis of errors

Last but not least, we are going to explore possible differences with regards to mean accuracies using a Bayesian ANOVA.

# create accuracies dataframe
df.acc = df.tsk %>% filter(expected != "neutral") %>% droplevels() %>%
  group_by(subID, diagnosis, phase, expected, difficulty) %>%
  summarise(acc = mean(acc)*100)

# check normal distribution
kable(df.tsk %>% group_by(subID, diagnosis, phase, expected, difficulty) %>%
  summarise(acc = mean(acc)*100) %>% 
  group_by(diagnosis, phase, expected, difficulty) %>%
  shapiro_test(acc) %>%
    mutate(
      sig = if_else(p < 0.05, "*", "")
    ))
diagnosis phase expected difficulty variable statistic p sig
ADHD prevolatile expected easy acc 0.7530551 0.0001005 *
ADHD prevolatile expected medium acc 0.7069895 0.0000240 *
ADHD prevolatile expected difficult acc 0.8478059 0.0031053 *
ADHD prevolatile unexpected easy acc 0.7684096 0.0001668 *
ADHD prevolatile unexpected medium acc 0.6030083 0.0000014 *
ADHD prevolatile unexpected difficult acc 0.6419123 0.0000039 *
ADHD volatile expected easy acc 0.7936914 0.0003993 *
ADHD volatile expected medium acc 0.8355305 0.0019035 *
ADHD volatile expected difficult acc 0.8366283 0.0019874 *
ADHD volatile neutral easy acc 0.7837217 0.0002813 *
ADHD volatile neutral medium acc 0.8443747 0.0027041 *
ADHD volatile neutral difficult acc 0.8929161 0.0215045 *
ADHD volatile unexpected easy acc 0.7557361 0.0001096 *
ADHD volatile unexpected medium acc 0.8294574 0.0015023 *
ADHD volatile unexpected difficult acc 0.8536698 0.0039446 *
ADHD postvolatile expected easy acc 0.8440215 0.0026661 *
ADHD postvolatile expected medium acc 0.9104130 0.0482689 *
ADHD postvolatile expected difficult acc 0.8661620 0.0066462 *
ADHD postvolatile unexpected easy acc 0.5959201 0.0000012 *
ADHD postvolatile unexpected medium acc 0.5555850 0.0000005 *
ADHD postvolatile unexpected difficult acc 0.7380254 0.0000621 *
BOTH prevolatile expected easy acc 0.7693266 0.0001720 *
BOTH prevolatile expected medium acc 0.8136402 0.0008245 *
BOTH prevolatile expected difficult acc 0.8606902 0.0052777 *
BOTH prevolatile unexpected easy acc 0.6126826 0.0000018 *
BOTH prevolatile unexpected medium acc 0.5222700 0.0000002 *
BOTH prevolatile unexpected difficult acc 0.8256146 0.0012958 *
BOTH volatile expected easy acc 0.8657282 0.0065251 *
BOTH volatile expected medium acc 0.7762210 0.0002173 *
BOTH volatile expected difficult acc 0.9051781 0.0377818 *
BOTH volatile neutral easy acc 0.7245658 0.0000409 *
BOTH volatile neutral medium acc 0.6932121 0.0000161 *
BOTH volatile neutral difficult acc 0.9250169 0.0966207
BOTH volatile unexpected easy acc 0.6507543 0.0000049 *
BOTH volatile unexpected medium acc 0.7485284 0.0000868 *
BOTH volatile unexpected difficult acc 0.7445783 0.0000765 *
BOTH postvolatile expected easy acc 0.7509063 0.0000937 *
BOTH postvolatile expected medium acc 0.8264958 0.0013403 *
BOTH postvolatile expected difficult acc 0.9085878 0.0443054 *
BOTH postvolatile unexpected easy acc 0.4955981 0.0000001 *
BOTH postvolatile unexpected medium acc 0.5505637 0.0000004 *
BOTH postvolatile unexpected difficult acc 0.7628110 0.0001384 *
COMP prevolatile expected easy acc 0.6048243 0.0000015 *
COMP prevolatile expected medium acc 0.6717248 0.0000087 *
COMP prevolatile expected difficult acc 0.9220654 0.0838957
COMP prevolatile unexpected easy acc 0.5904384 0.0000010 *
COMP prevolatile unexpected medium acc 0.5555850 0.0000005 *
COMP prevolatile unexpected difficult acc 0.8469958 0.0030052 *
COMP volatile expected easy acc 0.8053882 0.0006083 *
COMP volatile expected medium acc 0.7784453 0.0002345 *
COMP volatile expected difficult acc 0.8934097 0.0219914 *
COMP volatile neutral easy acc 0.6048243 0.0000015 *
COMP volatile neutral medium acc 0.7416438 0.0000696 *
COMP volatile neutral difficult acc 0.8721173 0.0085735 *
COMP volatile unexpected easy acc 0.6046337 0.0000015 *
COMP volatile unexpected medium acc 0.6449235 0.0000042 *
COMP volatile unexpected difficult acc 0.7733926 0.0001974 *
COMP postvolatile expected easy acc 0.6488874 0.0000047 *
COMP postvolatile expected medium acc 0.7015343 0.0000205 *
COMP postvolatile expected difficult acc 0.8164406 0.0009154 *
COMP postvolatile unexpected easy acc 0.4735522 0.0000001 *
COMP postvolatile unexpected medium acc 0.4225841 0.0000000 *
COMP postvolatile unexpected difficult acc 0.7202060 0.0000358 *
# rank transform the data
df.acc = df.acc %>%
  ungroup() %>%
  mutate(
    racc   = rank(acc)
    )

# run the ANOVA
if (!file.exists(file.path(brms_dir, "aov_acc.rds"))){
  aov.acc   = anovaBF(racc   ~ diagnosis * phase * expected * difficulty, data = df.acc)
  saveRDS(aov.acc, file = file.path(brms_dir, "aov_acc.rds"))
} else {
  aov.acc = readRDS(file.path(brms_dir, "aov_acc.rds"))
}

bf.acc = extractBF(aov.acc, logbf = T)
kable(head(bf.acc %>% arrange(desc(bf))))
bf error time code
diagnosis + expected + difficulty 46.41268 0.0176216 Fri Sep 5 15:14:41 2025 23cdf312fa7
diagnosis + difficulty 46.24773 0.0220083 Fri Sep 5 15:14:41 2025 23cd20502937
diagnosis + expected + difficulty + diagnosis:difficulty 45.01890 0.0189216 Fri Sep 5 15:14:44 2025 23cd306bedbc
diagnosis + difficulty + diagnosis:difficulty 44.85054 0.0153529 Fri Sep 5 15:14:43 2025 23cd186e77cb
diagnosis + phase + expected + phase:expected + difficulty 44.76534 0.1771458 Fri Sep 5 15:14:42 2025 23cd2da83298
expected + difficulty 44.41762 0.0165907 Fri Sep 5 15:14:41 2025 23cd4f1e8132
# print overall accuracy rates for all the effects included in the best model
df.acc %>% 
  group_by(diagnosis, difficulty, expected) %>% 
  summarise(mean_accuracy = mean(acc, na.rm = T), 
            sd_accuracy = sd(acc, na.rm = T))
## # A tibble: 18 × 5
## # Groups:   diagnosis, difficulty [9]
##    diagnosis difficulty expected   mean_accuracy sd_accuracy
##    <fct>     <fct>      <fct>              <dbl>       <dbl>
##  1 ADHD      easy       expected            94.1        6.87
##  2 ADHD      easy       unexpected          88.4       17.0 
##  3 ADHD      medium     expected            92.6        8.14
##  4 ADHD      medium     unexpected          90.7       14.3 
##  5 ADHD      difficult  expected            88.6        9.62
##  6 ADHD      difficult  unexpected          86.6       16.0 
##  7 BOTH      easy       expected            94.9        6.58
##  8 BOTH      easy       unexpected          92.4       12.8 
##  9 BOTH      medium     expected            93.9        7.07
## 10 BOTH      medium     unexpected          91.9       13.8 
## 11 BOTH      difficult  expected            87.9        9.95
## 12 BOTH      difficult  unexpected          83.0       18.6 
## 13 COMP      easy       expected            97.9        3.23
## 14 COMP      easy       unexpected          93.9       10.8 
## 15 COMP      medium     expected            96.3        5.34
## 16 COMP      medium     unexpected          94.7       11.4 
## 17 COMP      difficult  expected            89.8        8.68
## 18 COMP      difficult  unexpected          84.8       16.6
df.acc %>% 
  group_by(diagnosis) %>% 
  summarise(mean_accuracy = mean(acc, na.rm = T), 
            sd_accuracy = sd(acc, na.rm = T))
## # A tibble: 3 × 3
##   diagnosis mean_accuracy sd_accuracy
##   <fct>             <dbl>       <dbl>
## 1 ADHD               90.2        12.8
## 2 BOTH               90.7        12.8
## 3 COMP               92.9        11.1
df.acc %>% 
  group_by(expected) %>% 
  summarise(mean_accuracy = mean(acc, na.rm = T), 
            sd_accuracy = sd(acc, na.rm = T))
## # A tibble: 2 × 3
##   expected   mean_accuracy sd_accuracy
##   <fct>              <dbl>       <dbl>
## 1 expected            92.9        8.18
## 2 unexpected          89.6       15.2
df.acc %>% 
  group_by(difficulty) %>% 
  summarise(mean_accuracy = mean(acc, na.rm = T), 
            sd_accuracy = sd(acc, na.rm = T))
## # A tibble: 3 × 3
##   difficulty mean_accuracy sd_accuracy
##   <fct>              <dbl>       <dbl>
## 1 easy                93.6        10.9
## 2 medium              93.4        10.7
## 3 difficult           86.8        13.9

Accuracies were generally high, with a grand average of 91.25% accurate responses across diagnostic groups. Accuracies seems to have differed between diagnostic groups. Let’s do some plotting to figure out where our differences lie.

4.1 Plots

# rain cloud plot including all factors

df.acc %>%
  ggplot(aes(expected, acc, fill = diagnosis, colour = diagnosis)) + #
  geom_rain(rain.side = 'r',
boxplot.args = list(color = "black", outlier.shape = NA, show.legend = FALSE, alpha = .8),
violin.args = list(color = "black", outlier.shape = NA, alpha = .8),
boxplot.args.pos = list(
  position = ggpp::position_dodgenudge(x = 0, width = 0.3), width = 0.3
),
point.args = list(show_guide = FALSE, alpha = .5),
violin.args.pos = list(
  width = 0.6, position = position_nudge(x = 0.16)),
point.args.pos = list(position = ggpp::position_dodgenudge(x = -0.25, width = 0.1))) +
  scale_fill_manual(values = col.grp) +
  scale_color_manual(values = col.grp) +
  labs(title = "Accuracies per subject", x = "", y = "accuracy (%)") +
  facet_wrap(phase ~ difficulty) +
  theme_bw() + 
  theme(legend.position = "bottom", plot.title = element_text(hjust = 0.5), 
        legend.direction = "horizontal", text = element_text(size = 15))

# plot as line plot without the phases
df.acc %>%
  group_by(subID, diagnosis, expected, difficulty) %>%
  summarise(
    acc = mean(acc, na.rm = T)
  ) %>% 
  group_by(diagnosis, expected, difficulty) %>%
  summarise(
    acc.mn = mean(acc),
    acc.se = sd(acc)
  ) %>%
  ggplot(aes(y = acc.mn, x = difficulty, group = expected, colour = expected)) +
  geom_line(position = position_dodge(0.4), linewidth = 1) +
  geom_errorbar(aes(ymin = acc.mn - acc.se, 
                    ymax = acc.mn + acc.se), linewidth = 1, width = 0.5, 
                position = position_dodge(0.4)) + 
  geom_point(position = position_dodge(0.4), size = 5) + 
  facet_wrap(. ~ diagnosis) +
  #ylim(75,100) +
  labs (x = "difficulty", y = "accuracy") + 
  scale_fill_manual(values = col.cnd) +
  scale_color_manual(values = col.cnd) +
  theme_bw() + 
  theme(legend.position = "bottom", plot.title = element_text(hjust = 0.5))